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1 Partial differential equations and numerical methods 



1.1 Diffusion, convection and reaction 

We shall use the following notations for the so-called gradient and divergence operators, 
grad f(x,y,z) = I ^f(x,y,z 



' fi(x,y,z) 
f2(x,y,z) 
,h(x,y,z) 



o^f(x,y,z) 

ox ay az 



Q Q 

div I f 2 (x,y,z) I = — + — f 2 (x,y,z) + — f 3 (x,y,z) 



We consider a chemical species in a three-dimensional region £1, and denote its concentration 
(in mass per unit volume), at the time t and location p = (x, y, z) in 17, by u(p, t) = u(x, y, z, t). 

We denote the rate of generation of the chemical species per unit volume, due to some 
chemical reaction, by r(p,t,u(p,t)). This rate thus depends on the position p, time t and actual 
concentration u(p, t). The rate can have a positive or negative value, and be zero for a nonreacting 
species. If there would be no flow of the species through f2, e.g. if 0, consists of a substance 
impermeable to the species, the concentration would satisfy the differential equation 

d 

— u(p,t) = r(p,t,u(p,t)). 

But, in many cases of practical interest (fluids or gasses) we have to consider a flux through 
f2 influencing u(p,t). We represent this flux by a vector F(p,t) £ M 3 with components Fi(p,t), 
F2(p,t), F3(p,t). This vector has a direction coinciding with the direction of the flux, and its 
Euclidean length equals the amount of species flowing per time unit through a unit surface that is 
perpendicular to the flux (at the position p and time t). In order to determine the effect of the flux 
F on -§iu(p, t) we consider a small cube AO in with center at p = (x, y, z). Let AQ be bounded 
by surfaces perpendicular to the x—,y— and z— axes with distances from p equal to 4p, ^n, 
respectively. By simple geometric arguments it is possible to determine the effect of the flux on the 
increase of the amount of species, per time unit, in the cube. This effect is approximately equal to 



Ax , Ax 



2 ' ' ^ 2 



Fi(x- — ,y,z,t) -F 1 (x+ — ,y,z,t) AyAz + F 2 (x,y - —,z,t) - F 2 (x,y + —,z,t) 



+ 



that is to 



Az . ^ . Az . 

F3(x,y,z- —,t) - F 3 (x,y,z + ~Y,t) 



-div F(p,t) ■ Ax AyAz. 



2 ' ' ' * v 2 
AxAy, 



AxA; 



Therefore the effect of the flux on the increase of u(p,t), per time unit, equals —div F(p,t). The 
total effect of reaction and flux is expressed in 

d 

— u(p, t) = -div F(p, t) + r(p, t, u(p, t)), 



which we shall call the conservation of mass equation. 
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In general the flux vector F(p, t) can be decomposed, 

F(p,t) = u(p,t) -V(p,t) + W(p,t). 

Here the first term in the right-hand member is the convection flux of the species. It is caused 
by a general flow field V(p, t) G M 3 which is independent of the concentration of the chemical 
species. The amount of the species transported with the motion V(p, t) is represented by the 
vector u(p,t) ■ V(p,t). The second term is called the diffusion flux. It is due to the molecular, 
thermal agitation, and is also present in fluids or gasses at rest. According to Ficks' law 

W(p,t) = -k(p,t) ■ grad u(p,t), 

where k(p, t) is called the diffusion coefficient at p and at time t. 

Substituting the above expressions in the conservation of mass equation we obtain the partial 
differential equation 

d 

(1.1.1) ~Q t U ( p i *) = dlV ™ *) grad *)] ~ dlV ™ *) V ( p i *)] + *' U (P' *))■ 

It is called a transport equation, or also a diffusion- convection-reaction equation, for u(p,t). 

In practical situations one may be interested in computing the solution u{p,i) to (1.1.1) for 
p G fl, t > 0, under an initial condition of the form 

(1.1.2) «(p,0) = /(p) for p en, 

with given function /(p). In general, the relations (1.1.1), (1.1.2) are supplemented by conditions 
for t > at the points p belonging to the boundary <9f2 of 0. These boundary conditions are usually 
of one of the following three types, 

(1.1.3a) u{p,t) = g{p,t), 

(1.1.36) —u{p,t)=g{p,t), 

d 

(1.1.3c) —u(p,t) = c{p,t) ■ \g{p,t) -u(p,t)\. 

Here g(p, t) and c(p, t) are given functions. Further, n denotes the outward normal to dVL at p, with 
length 1, so that -J^-u(p,t) is equal to the inner product of n and grad u(p,t). 

Condition (1.1.3a) is called a Dirichlet boundary condition. It occurs if there is free contact 
between the interior and exterior of f2 at p £ dQ, and the concentration outside 0, is known to be 
equal to uq (p, t) . If the medium in which the chemical species is solved is the same — outside and 
inside of 30 — we have g(p, t) = uo(p, t). 

Condition (1.1.3b) is called a Von Neumann boundary condition. E.g. if the boundary <9f2 
is impermeable to the chemical species, and the normal component V n (p,t) of V{p,t) at p 6 d£l 
vanishes, we have (1.1.3b) with g(p,t) = 0. 

Condition (1.1.3c) is a mixed condition. E.g. if the boundary is semipermeable with regard 
to the chemical species, we may have c(p,t) > and g(p,t) as in (1.1.3a). 

In cases where -^u(p,t) = 0, or -^u{p,t) = -^-u(p,t) = 0, throughout f2 and for all t > 0, it 
may be convenient to consider the concentration as a function only of the variable (x, y, t) or (x, t), 
respectively. In this way one can arrive e.g. at the following four lower-dimensional versions of the 
above. 
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(1.1.4) —u{x,t) = —u{x,t), u(0,t) = g (t), u(M)=0i(t), 
u (x, 0) = /(x), where < x < 1, t > 0. 

This is the pure diffusion equation with Dirichlet boundary conditions in 1 dimension. 

(1.1.5) ^u{x,t) = -^[b(x)u(x,t)]+c{x,t), u(0,t) = 0, 

n (x, 0) = /(a:), where < x < 1, t > 0. 

This is a one-dimensional convection-reaction equation, where the rate c(x, i) of the chemical reac- 
tion is independent of the concentration u(x,t). 



(1.1.6) —u(x,t) = — [a(x)—u(x,t)] + — [b(x)u(x,t)] + c(x)u(x,t) + d(x), 

d 

u(0, i) = gi(t), fa; u 0-, t) = °> °) = /<>)> where < x < 1, t > 0. 

This is a one-dimensional diffusion-convection-reaction equation. In case 6(0) < 0, the convection at 
x = is from left to right, and the boundary point x = lies upwind, sometimes called upstream, of 
the calculation domain. Correspondingly, the condition u(0, t) = g(t) is called an inflow boundary 
condition. In case 6(0) > 0, the same condition would be called an outflow boundary condition. 

In the following example we consider the square Q = {(x, y) : < x < 1, < y < 1}, with 
boundary dQ = {(x, y) : (x, y) G Q, and x(x — 1) y(y — 1) = 0}. 



d d 2 d 2 

(1.1.7) di U ( X > y ' = 'dx 2 ^ 30 ' y ' ^ + dy 2 ^' V: ^ f ° r ^' G °' 

u{x,y,t) = g(x,y,t) for (x,y) e dQ, 

u{x, y, 0) = f(x, y) for (x, y) G fi, 

where t > 0. 

This is a pure diffusion equation with Dirichlet boundary conditions in 2 dimensions. 



1.2 Semi-discretization by finite difference methods 
1.2.1 Basic finite difference approximations 

The finite difference method can be used for obtaining numerical approximations to the solutions 
of the partial differential equations reviewed above under given initial-boundary conditions. The 
method essentially consists in replacing the partial derivatives in the original problem by finite 
difference quotients so as to obtain a related problem which is easier to solve than the original 
one. In this section we focus on such replacements for the derivatives with respect to the space 
variables (x, y and z) only. This process, where derivatives with respect to t are not affected, is 
called semi-discretization. 

We start by listing various difference approximations to the derivatives of a function v(x) 
defined for < x < 1. We assume v(x) to have continuous derivatives of the orders occurring in 
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the following expressions, and for integer values p > 1 we denote by v^(x) the p-th derivative of 
v(x). Further, we use the notation 



l-y(p) I = max \y( p \a 

1 lo ° 0<x<l' 



)|- 



All arguments of v, occurring below, belong to the interval [0, 1], and 5 = Ax denotes a positive 
increment of the variable x. The following six relations are valid: 

(1.2.1a) v'(x) = <T>(:e + 5)- v(x)] + R, 

with \R\ < -{v^loo (forward difference approximation). 

(1.2.16) v'{x) = 5' 1 [v(x) - v(x - 5)] + R, 

with \R\ < -{v^loo (backward difference approximation). 

(1.2.1c) v'{x) = {25)- 1 [v(x + 5)- v(x - 6)] + R, 

5 2 

with \R\ < — \v^\oo (centered approximation). 
(1.2.2a) v"(x) = 5- 2 [v{x-5) -2v(x) +v(x + S)] + R, with \R\ < — {v^l^. 

(1.2.26) v"{x) = ^6~ 2 [-v{x - 25) + 16u(x -8)- mv(x) + 16v(x + 5)- v(x + 25)} + 

90 1 

(1.2.2c) j^[ v "( x ~ s ) + 10u"(x) + v"(x + 5)} = 5- 2 [v(x -5)- 2v(x) + v(x + 5)} + R, 

5 A 

with \R\ < ^^l w ^ 6 ' ) |oo (Numerov's approximation). 

The above upper bounds for \R\ can be obtained by expanding R in powers of 5, with 
the use of Taylor's formula applied to the function v. In order to arrive at the bounds in (1.2.2b), 
(1.2.2c) one has to express the remainder term associated with the Taylor polynomial as an integral 
involving the 6-th derivative of v. Further, in proving (1.2.2c), Taylor's formula should as well be 
applied to the function v" . 

The left-hand members in the formulas (1.2.1), (1.2.2) can be approximated by the corre- 
sponding right-hand members in which the residual R is suppressed. 



+ R, with \R\ < —\v { , . 



1.2.2 Finite difference approximations using non-equidistant arguments 

It is also possible to establish finite difference approximations with nonuniformly distributed argu- 
ments. For instance, let Oi > 0, 52 > and consider the expression 

P ■ v(x - Si) + a ■ v(x) + 7 • v(x + 5 2 ), 



where a, (3, 7 are coefficients which are still to be determined. By Taylor's formula this expression 
equals 



{a + (3 + j)v{x) + ( 7 5 2 - pS x )v'{x) + j(5 2 ) 2 + 



v"(x) 



R, 



where 
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\R\< h\(S2) 3 + \P\(5i 



1 

J 6 



«'"|oo. 



Suppose one wants to approximate v"(x). Then the requirements on a, (3, 7 should be that 
a + /3 + 7 = 0, 7*2-^=0, 7(5 2 ) 2 + /?(5i) 2 = 2, 



which yields 



a 



<5i<!>2 ' 

From the above it can be seen that 



5i(S 1 +6 2 y 
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(1.2.3) 



2 2 2 
i>"(x) = -—- ^-v(x — Si) — — — v(x) + -—- ^-v(x + do) + R, 

with < - ■ max(<5i, <5 2 ) ■ \v"'\oo- 



1.2.3 Semi-discretization in one space variable 

In the process of semi-discretization by finite differences we choose gridpoints p in the calculation 
domain, and we approximate the values u(p, t) by quantities U(p, t) for all gridpoints p. 

For instance, consider the approximation of the solution to (1.1.4) on a so-called nonuniform 
grid using formula (1.2.3). We choose 5\ > 0, 82 > 0, . . . , <5 s+ i > with 8\ + S 2 + ■ ■ ■ + 5 s +i = 1) 



and define gridpoints p 



by x = 0, x j+ i 



+ Sj + i (for j = 0, 1, . . . , s) . In view of (1.2.3) 



we define approximations Uj(t) = U(xj,t) ps u(xj,t) by requiring, for 1 < j < s, 



<^(5j +<5j-+i) 



^■(t) + 



S j+1 (Sj +5 j+1 ) 



U j+1 (t), U j (0) = f(x 



and by putting Uo(t) = go(t), U s+ i(t) = g\{t). These requirements can be cast in the compact 
form 



(1.2.4) 



d_ 

dt 



U(t) = AU(t) + r(t) for t > 0, 17(0) = u . 



Here A is a tridiagonal s x s matrix, r(t) and tio are given vectors in the s-dimensional complex 
vector space C s , and U (t) £ C s has components Ui(t), Uz(t), ■ ■ ■ , U s (t). The solution to the system 
of ordinary differential equations in (1.2.4) thus provides an approximation to the solution of the 
partial differential equation in (1.1.4). 

The use of variable increments 5j can be advantageous in cases where the smoothness of the 
true solution u(x,t) varies significantly when x runs through [0,1]. In such cases one may choose 
Sj especially small in those parts of [0, 1] where u(x, t) is expected to be nonsmooth. 

If the smoothness of u(x, t) is not expected to vary strongly with x it may be recommended 
to use a so-called uniform grid, i.e. Xj+i = Xj + 5 with 5 = (s + 1) _1 for j = 0, 1, . . . , s. Similarly 
as in the case of a nonuniform grid, we arrive at a semi-discrete problem of the form (1.2.4). But 
now we can use simply (1.2.2a) instead of the more general formula (1.2.3). In this case the matrix 
A is of the simple form 

-2 1 
1 -2 I 

A 
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More accurate approximations to the solution of (1.1.4) may be obtained, on a uniform grid, 
if one replaces the second order derivatives -^^u(xj,t) by the finite difference expression in (1.2.2b) 
for j = 2, 3, . . . , s — 1 and only by the expression from (1.2.2a) for j = 1, s. In this way one again 
arrives at a semi-discrete problem of the form (1.2.4) but now with a pentadiagonal matrix A. 

Also Numerov's approximation (1.2.2c) can be used to construct a semi-discrete version of 
(1.1.4). On a uniform grid, where Xj = j • 5, 5 = (s + 1) _1 , one arrives, for 1 < j < s, at the 
requirements 

^p/;_i(t) + 1017<(t) + U' j+1 (t)] = S-^Uj-^t) - 2Uj{t) + U j+1 (t)], Uj(0) = f( Xj ). 
By using Uo{t) = go(t), U s+ i(t) = gi{t) these requirements can be formulated in the compact form 
A 1 jU{t) = A U{t) + q{t) fort>0, U(0) = u o , 

with known vectors q(t), uq € C s , and tridiagonal s x s matrices Ao,Ai. Multiplying by (Ai) _1 
we see that this initial value problem for U(t) is of the form (1.2.4) with A = (A 1 )~ 1 A , r(t) = 
(A 1 )-'q(t). 



1.2.4 Semi-discretization in more than one space variable 

In the case of two or three space variables we can proceed similarly as above. 

As an illustration we consider the semi-discretization of problem (1.1.7) using the finite 
difference approximation (1.2.2a). Choose an integer N > 1, 5 = (N + 1) , and define the 
gridpoints pj k = (j5, k5) for j = 0, 1, . . . , TV + 1 and k = 0, 1, . . . , N + 1. Replacing the partial 

derivatives J^u(pjk, t), -§^u(pjk, t) occurring in (1.1.7) (when x = jS, y = k5) by the finite 
difference approximations 

$~ 2 [u(Pj-i,k,t) - 2u(p jk ,t) + u{p j+1;k ,t)], 5~ 2 [u{p j;k -i,t) - 2u(p jk ,t) + u(p jtk+1 ,t)], 

respectively, one obtains the system of ordinary differential equations 

( d 

—U(p jk ,t) = 5~ 2 [U{pj- ltk ,t) + U(p j+li k,t) + U(p jtk - 1 ,t) + U(pj, k+1 ,t) - 4U(p jk ,t)], 
where j = 1, 2, . . . , N and k = 1, 2, . . . , N. 

In the right-hand members of this system we set U (p mn , t) = g(p mn ,t) for all p mn 6 dQ (cf. (1.1.7)). 
Moreover, we supplement the system by the initial conditions 

U(j)jk, 0) = ffyjk), where j = 1, 2, . . . , N and k = 1, 2, . . . , N, 

in conformity with the initial condition in (1.1.7). 

For t > we introduce a vector U(t) G C s , with s = N 2 , the components of which are, 
in some fixed order, the values U(pj k ,t) (with 1 < j < N, 1 < k < N). The above system of 
A^ 2 ordinary differential equations, together with the corresponding initial conditions, can again be 
cast in the form (1.2.4). The components of uq are equal to values f(pjk), an d the components 
of r(t) depend on values g(p mn ,t)- Further, each row of the matrix A contains at most 5 nonzero 
entries. The entries on the main diagonal equal — 4/(5 2 , and the nonzero off-diagonal elements are 
equal to 1/5 2 . 
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1.2.5 Semi-discretization of diffusion and convection phenomena 

1.2.5.1 Approximating diffusion and convection terms 

We now consider in some detail various finite difference approximations to the diffusion and con- 
vection terms occurring in (1.1.6). 

We shall use, for any real £, the notations 



£ if e > 0, 
otherwise 



and 



£ if € < 0, 
otherwise 



One easily sees that 

e = e + +r, iei=e + -r, (-o~ 



■(0 _ , (-0" 



■(0- 



In the following we consider an increment 5 > and given function a(x), b(x), v(x). Further, 
we denote by e{x) certain parameter values to be specified below. 

We consider the following four finite difference approximations. 

(1.2.5a) (b(x)v(x))' = {/3 ■ v(x - 5) + a ■ v(x) + 7 • v(x + 5)} + R, with 

7 = 7 ( x ) = 5" 1 [;U(x + 5/2) + e(s + 5/2) + 5/2)] + - ^6(x + 5/2))] , 

= p{x) = 5' 1 --b(x - 5/2) + e(x - 5/2) (-6(x - 5/2) - [b(x - 5/2)]") 
a = a {x) = -(3(x) - 7 (x) + 5 _1 [-6(x - 5/2) + b(x + 5/2)}. 



(1.2.56) (b(x)v(x))' = {/? ■ v(x - 5) + a • v(x) + 7 • v(x + 5)} + R, with 

7 = 7 ( x ) = J" 1 [^6(x + 5) + e(x + 5) ([6(x + 5)]+ - ^6(x + 5)) 

= /?( x ) = 5' 1 --b(x -5)+ e(x - 5) (-6(x -5)- [b(x - 5)]") 
a = a(x) = 5" 1 £(x)([6(x)]- - [b(x)] + ). 



(1.2.6a) 



(a(x)v'(x))' = {(3 ■ v(x - 5) + a • v{x) + 7 • v(x + 5)} + i?, with 
7 = 7 ( x ) = 5~ 2 a(x + 5/2), 
P = P(x) = 5" 2 a(x-5/2), 
a = a(x) = —(3(x) — 7 (x). 



(1.2.66) (a(x)u' (x))' = {/9 • v(x - 5) + a ■ v(x) + 7 • u(x + 5)} + with 

7 = 7 ( x ) = ^5" 2 [a(x) + a(x + 5)], 

/3 = /3(x) = i5- 2 [a(x-5)+a(x)], 
a = a(x) = —(3(x) — 7 (x). 

For a(x), t>(x) sufficiently smooth the residuals R in (1.2.6a), (1.2.6b) can be seen to be 0(5 2 ) - 
which is similar to the residual in (1.2.2a). If e(x) = and 6(x), v(x) are sufficiently smooth, then 
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also the residuals R in (1.2.5a), (1.2.5b) can be seen to be 0(5 2 ) — which is similar to the residual 
in (1.2.1c). But, in general, if the e-values in (1.2.5a) or (1.2.5b) are nonvanishing parameters, then 
the corresponding residuals are only 0(5). 

In order to explain why parameter values different from zero may still be useful we focus on 
the physical meaning of the convection term in (1.1.6). In fact, -^[b(x)u(x,t)] can be interpreted 
as the rate at which the concentration of a chemical species increases, at the position x and time 
t, due to the 1-dimensional flow field V(x) = —b(x) (cf. Section 1.1). 

For the time being assume b{x) > 0. Then the flow is from right to left, and therefore the 
values u(x,t') with t' > t will be influenced by the so-called upwind values u(x',t) with x' > x. If 
one wants to use a finite difference approximation to ^[b(x)u(x , t)] for determining u(x,t') with 
t' > t, it thus seems quite natural to choose the forward difference approximation (1.2.1a) (with 
v replaced by b ■ u) and not (1.2.1b) or (1.2.1c). The last two approximations would force u(x,t') 
with t' > t to depend on concentration values u(x',t), with x' < x. But, these values should have 
no influence on u(x,t'). Similarly, for b(x) < it is natural to use (1.2.1b) (with v replaced by 
b-u). 

Putting e(x — 5) = e(x) = e(x + 5) = 1 in (1.2.5b), this formula reduces to (1.2.1a) (with v 
replaced by b- v) if b(x — 5) > 0, b(x) > 0, b(x + 5) > 0; and it reduces to (1.2.1b) (with v replaced 
by b • v) if b(x — 5) < 0, b(x) < 0, b(x + 5) < 0. With all parameter values equal to 1 formula 
(1.2.5b) is called an upwind finite difference approximation. Although, for small S > 0, its accuracy 
is lower than of (1.2.5b) with all e-values equal to 0, it is the more natural approximation from a 
physical point of view. In practical applications the upwind parameters e(x) are often chosen such 
that < e(x) < 1. 

Arguments similar to the above apply to the finite difference approximation in (1.2.5a). 
Again the values e(x — 5/2) and e(x + 5/2) are called upwind parameters. 

1.2.5.2 Semi-discretization of problem (1.1.6) 

We shall construct a semi-discrete version of (1.1.6). We choose an integer s > 2, we put 5 = 
x\ = A • 5, and we define a\ = a(x\), b\ = b(x\), c\ = c(x\), d\ = d(x\). Approximations 
Uj(t) w u(xj,t) (for j = 1,2, . . . ,s) can be obtained by solving 

(1.2.7a) J t U ^ = AU ^ + ^ f ° r * ~ °' ^ = U °' 

Here U(t) G C s has the components Uj(t) (for j = 1, 2, . . . , s), the initial vector u G C s has the 
components f(xj) (for j = 1,2,... , s), and the vectors r(t) G C s depend on dj, g(t). The s x s 
matrix A is of the form 

(1.2.76) A = A 2 + A l +A , 

where A 2 , A\, Aq correspond to the first three terms in the right-hand member of the partial 
differential equation in (1.1.6). 

We define Aq to be the diagonal matrix 

(1.2.8) A = (ajk) with ctjj = Cj, ctjk = (for 1 < j < s, 1 < k < s, j ^ k). 

The s x s matrices A\, A 2 are both tridiagonal; they can be constructed by using (1.2.5a), 
(1.2.6a) or (1.2.5b), (1.2.6b). In all cases we denote the entries of Ai, A 2 by ctjk (for 1 < j < s, 1 < 
k < s), and we define 

a l,j+i = l( x i) (! < 3 < s ~ 1), 

a 3,j-i = P( x l) ( 2 < i < s - 1), a SjS _i = 0(x s ) + 7(x s ), 
a l,j = a ( x l) (! < 3 < s )- 
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Here a(x), f3(x), 7(2) are as in (1.2.5a), (1.2.6a), (1.2.5b) or (1.2.6b). In the above expression for 
a s s _i we have included the term j(x s ) so as to express the boundary condition -^-u(l,t) = (see 
(1.1.6)). 

Basing our construction on (1.2.5a), (1.2.6a) we arrive at the following definitions (1.2.9a), 
(1.2.10a), and basing it on (1.2.5b), (1.2.6b) at the definitions (1.2.9b), (1.2.10b). 



(1.2.9a) 



M = (ajk) with , for j = 1, 2, . . . , , 



5- 1 



L+e j+ ilb+ + >--b 



1 



(1.2.96) 



M = (a jk ) with, for j = 1, 2, . . . , s, 



a j,j+i 



= ,5- 1 



a i>j-i = 5 1 



-\bj-i+ £ j-i(\ b j-i ~ b j-i 



a 



3,3 



S-'sjlbJ-bt}. 



(1.2.10a) 



A 2 = (a jk ) with 

a jJ+1 = S~ 2 a j+ i ( for 1 < j < s - 1), 



a 



3:3 



a 



3,3 



-i = 5 2 [aj_i + 6 jt8 a a+ i] ( for 2 < j < s), 
= -(T 2 [a-i + a-+i] ( for 1 < j < s). 



(1.2.106) 



-42 = (cKjfc) with 



O: 



3,3+1 = j 5 2 [ a i + a i+i] ( for 1 < j < s — 1), 

= -<5~ 2 [a.,_i + a, + <5j, s (a a + a s+ i)] ( for 2 < j < s), 



a 



o a j-i + «j + - fl 



( for 1 < j < s). 



In the above four definitions we have used the Kronecker index 5jk defined by 



Sjk = 1 for j = k, and 5jk = for j / k. 



Further, the values a\ = a(x\), b\ = b(x\) which occur above with x\ > 1 should be interpreted 
as obtained by a mathematical extension of the given functions a(x), b(x) defined for < x < 1. 
These values a\ , 6a need not to have a relation to any physical quantities outside the computation 
domain [0, 1]. In (1.2.9a) and (1.2.9b) the quantities E\ are upwind parameters, with < E\ < 1, 
and the quantities a SjS+ i, a^o have been introduced only for notational convenience. 
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1.3 Notes and remarks 

For further discussion of the physical processes considered in Section 1.1, see e.g. Crank (1975), 
Hirsch (1988), Patankar (1980) and Zlatev (1995). In these references the transport equation (1.1.1) 
is discussed as well as various boundary conditions like those in (1.3). 

There exists a vast literature on finite difference methods. The majority of the material 
on finite difference approximations in Section 1.2 can be found in one of the numerous books on 
this subject. We mention, in particular, the classical books Forsythe k, Wasow (1960), Richtmyer 
& Morton (1967) and the more recent works by Hirsch (1988), Meiss & Marcowitz (1981), Roos, 
Stynes & Tobiska (1996), Shashkov (1996), Strikwerda (1989), Thomas (1995), Thomee (1990). 

For the use of formula (1.2.2b), in deriving a semi-discrete version of a (nonlinear) diffusion- 
reaction problem with Dirichlet boundary conditions, see Stys & Stys (1991). 

Upwinding, as discussed in Section 1.2.5, is related to what some authors call the addition 
of artificial diffusion. For relevant literature, see e.g. Griffiths, Christie & Mitchell (1980), Hirsch 
(1988), Patankar (1980), Roos, Stynes & Tobiska (1996), Strikwerda (1989). 

Semi-discretization of (1.1.1) can be achieved by methods different from those described in 
Section 1.2, in particular by methods based on finite volumes, finite elements or (pseudo) spectral 
approximations. For finite volume methods see e.g. Hirsch (1988), Patankar (1980), Roos, Stynes & 
Tobiska (1996); for finite element methods see e.g. Hirsch (1988), Oden & Reddy (1976), Quarteroni 
& Valli (1994), Roos, Stynes & Tobiska (1996), Strang &; Fix (1973); and for pseudo spectral 
methods see e.g. Canuto, Hussaini, Quarteroni & Zang (1988), Fornberg (1996), Gottlieb & Orszag 
(1977), Quarteroni k Valli (1994). 
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2 Analysis 

2.1 Basic material 
2.1.1 Definitions 

Let (pit) denote a real function of a real variable t. This function is said to be isotone if <p{t\) < fit^) 
for all t\ < t 2 in the domain of definition of p. Conversely, if (p(ti) > tpfe) for all such ti,t2, then 
the function is called antitone. 

Let / denote any scalar function. If the A;-th derivative of / exists, it will be denoted by 
f k \ We define 

/ (0) = /• 

Let 7 = a + if3 = pe 10 denote any complex number, with real a, f3, 9 and p > 0. Then we 

write 

a = Re (7), (3 = Im (7), 7* = a — i/3, 

to denote the real part, imaginary part and the complex conjugate of 7, respectively. Further, if 
7 / 0, we define Arg(7), the principal value of the argument, to be equal to 6 with — ir < < ir. 
We define Arg (0) = 0, so that 

— 7r < Arg (7) < 7r for all 7 G C. 

Let 7 be any complex number, and p > 0. We denote the circ/e and rfis/c in the complex 
plane with center 7 and radius p by 

T[ 1 ,p\ = {C: |C-7|=P}, D[ 1 ,p] = {(: \( - 7 | < p}, 

respectively. 

Let V be an arbitrary subset of the complex plane. The intersection of all convex sets W 
with V C W C C is called the convex hull of and is denoted by 

conv V. 

Note that conv V is the smallest convex set in C containing the set V. The interior, closure and 
boundary of V are denoted by 

int V, cl V, dV, 
respectively. The distance from ( £ C to V is 

d(C,*0= mf{|C-r/| : 7/ G K}- 
The left half-plane is denoted by 

C_ = {C : (GC with Re(C) < 0}, 
and for 5 > the wedge W(S) is defined by 

= {C : C G C satisfies |Axg(— C)| < 6}. 
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Whenever m and n are integers with n < m, we adopt the conventions 

n n 

max • • • = 0, V ■ ■ ■ = 0, [T = 1, 

m<j<n * — ' 

j = rn j=m 

and we put 

0! = 1, 0° = 1. 

2.1.2 Theorems 

Let i?(C) = P(()/Q(() be a given rational function, where P((), Q(() are polynomials of degrees 
m and n, respectively. We assume that -P(C) and Q(() have no common zeros. 

Suppose Ci is a zero of Q(() with multiplicity k > 1. Then £i is called a poZe o/ order k of 
i?(C). We have Q(C) = (C - Ci) fe <5i(C), with Qi(C) of degree (n-/c) and Qi(Ci) / 0. By expanding 
P(C) and Qi(C) hi powers of (C — Ci) we see that, for C / Ci an d |C — Ci| sufficiently small, 

R(Q = a_ fc (C - (i)- k + ■ ■ ■ + a_i(C - Ci)" 1 + «o + «i(C - Ci) + a 2 (C - Ci) 2 + • • • ■ 

The expression a_fc(C — (i)~ k H + «-i(C — Ci) -1 is called the principal part of i?(C) at Ci, and 

a_i is the residue of i?(C) at Ci- 

By subtracting from R(() the principal parts corresponding to all different zeros of Q(C) we 
are left with a function f(() that is holomorphic on C. There are constants K, p such that for all 
( with \(\ > po we have the inequality 

i/(oi <K.\cr~ n . 

Further, for all ( £ C the value f(() can be represented as the sum of a convergent power series 

/(C) = 7o+7iC + 72C 2 + ---- 
Consequently, for any p > p , integration along the positively oriented circle T[0,p] yields 



/(C) 



<K-p 



m — n — k 



By letting p — > oo, we see that jk = for all k > m — n. We thus arrive at the following. 

Theorem 2.1.1 {Partial fraction decomposition) Let P(C), Q(C) be nonzero polynomials of de- 
grees m, n, respectively. Let Ci, C2> • • • > Cq be the zeros of Q((), with multiplicities fci, &2, . . . , k q , 
respectively. Then, for ( ^ (i, . . . , ( q , we have 

^=R (() + R 1 (() + ...+R q (0. 

Here i?o(C) 1S a nonzero polynomial of degree (m — n) if m > n, and i?o(C) = if m < n, so that 
we can write 

m — n 
1=0 

Further, for 1 < j < q, the quantities Rj(C) are ^ ne principal parts of P(()/Q(() at (j, so that we 
can write 



Rj(Q = X>-ij " Q~ l ( for 1 < j < q). 



1=1 
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Without proof we state the following important theorem concerning the factorial function 

n\ 

Theorem 2.1.2 (Stirlings formula) For each integer n > 1 there is a number 9, with < 9 < 1, 
such that 



n!=(^)V2^exp(— ). 



2.2 Estimates for the arc length of curves in the complex plane 
2.2.1 The arc length 

Let a < p. Suppose g(t) and h(t) are piecewise continuously differentiable real functions defined 
for a < t < (3. The function f(t) = g{t) + ih(t) defines a curve V in the complex plane. The length 
of the curve V can be specified by 



(2.2.1) L= I \f'(t)\dt. 

J a 



For each t G [a, (3] there is a real uj such that g'(t) = \f'{t)\ cosuj, h'(t) = \f'{t)\ sinw, which 

implies 



/ 1 27T / 1 27T 

/ \g'(t)cos9 + ti(t) sin 6\d6 = \ cos{u - 0)\ ■ \f{t)\d6 = 4|/'(t)|. 

Jo 

We thus arrive at the representation 

r-2-K r f3 



L=i^ |y \g' (t) cos6 + ti{t) sin 0\dt}( 



For each the quantity 



ft d 

L{6)= — {g (t) cos 9 + h(t) sin 6} 

J a 



dt 



is equal to the length of the projection of the curve T onto the straight line passing through the 
origin with angle 9 to the real axis. Therefore, we have proved the following expression of Cauchy 
for the length of a curve in terms of the length of its projections: 

(2.2.2) L = \ \ L(9)d6. 



4 
^0 



2.2.2 The image of a circle under a rational function 

We define a rational function R(C) = P(C)/Q(C) to be of order s if P(C), Q(C) are polynomials of 
degrees not exceeding s. We shall estimate the arc length of the image of a circle r[7, p] under a 
rational function of order s. This length can be written in the form 

l = r -^(7 / {Rhonda 

Jo dt Jt[i,p\ 
We shall relate the length L to the quantity 



M = max {|12(C) | = C lies on T[ 7 , p]}. 
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The main result of this subsection, Theorem 2.2.2, states that 

L < 2ns -M, 

whenever R(Q is a rational function of order s without poles on T[y, p\. 

This result is equivalent to the following remarkable fact. Denote the quantities L, M 
corresponding to the rational function P(() = (( — 7) s by Lp, Mp, respectively. Since Lp = 
2-ir s ■ p s , Mp = p s , we see that, for all R(() under consideration, the ratio L/M will never exceed 
the ratio Lp/Mp. 

Our proof of Theorem 2.2.2 will rely on the following lemma. 

Lemma 2.2.1 {On the arc length of the image of a circle under a rational function) Let R(() = 
P(() /(3(C) be a rational function of order s without poles on T [7, p] . Denote by M{9) the maximum 
value of the projection of R(() (where £ runs through T[y, p\) on the straight line passing through 
the origin with angle 6 to the real axis, i.e. 

M{6) = max {lie (e~ w R(()) : ( lies on r[7,/>]}. 

Then the length L of the image of T[y, p] under the mapping R satisfies 

r-2-IT 

L < s M{0)d6. 
Jo 

Proof. 1. With no loss of generality we assume 7 = 0, p = 1. In view of Cauchy's expression (2.2.2) 
it is sufficient to prove 

(2.2.3) L(9) <2s[M(6) + M(e + Tr)} for < 9 < 2tt, 

with 

fit) = g{t) + ih(t) = R(e H ), < t < 2vr. 
2. Let G [0, 2tt] be given, and write 

F(t) = g(t) cos 6 + h(t) sin 9. 

Since L{9) = J^\F'(t)\dt, we can assume with no loss of generality, that F'(t) does not vanish 
identically on [0, 2ir]. Therefore 

k t 

L(0) = El/' F '(*H' 

where to = < ti < . . . < tfc = 2tt and F'(t) / on each open interval (tj-i,tj). 
We introduce the notations 

a= min Fit), b= max F(t), 

0<t<27T 0<t<27T 

and we denote the range of values F(t) obtained when t runs through the interval (tj-i,tj) by 
F(tj-i,tj). Defining 

<Pj(x) = 1 for x G F(tj-i,tj) and <fj(x) = for x G [a, 6]\F(t J -_i, tj), 
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we have 

k r b r b k 



tpj{x)dx = / J2ipj(x)dx. 

7 = 1 



Below we shall prove 

k 

(2.2 A) ^ (pj(x) <2s for a < x < b. 

3 = 1 

Applying (2.2.4) we obtain 

L(0) < 2s(b -a) = 2s 



max Re(e" ie R(^)) - min Re(e- W R(e [t )) 



L 0<t<2n v ' 0<t<2n 

which proves (2.2.3). 

3. In order to prove (2.2.4) we assume x £ [a, b] to be given, and we write C = e 1 * for 
< t < 2it. A straightforward calculation shows that 

(2.2.5) F(t) = x 

is equivalent to 

e- ie P(o[Q(or +^[p(c)]*q(o - 2xQ(o[Q(or = o 

(where 7* denotes the complex conjugate of 7). Multiplying the latter equality by ( s we arrive, in 
view of C* = C _1 ) at a relation 

p(0 = 0, 

where p(£) is a polynomial of a degree not exceeding 2s. Moreover p(C) does not vanish identically 
(since F'(t) does not). Therefore, there exist at most 2s different values t £ (0, 2ir) with (2.2.5). 
This implies that x G F(tj-\,tj) for at most 2s different values j, so that (2.2.4) is fulfilled. □ 



Theorem 2.2.2 (On the arc length of the image of a circle under a rational function) Let R(C) 
be a rational function of order s without poles on F[y, p]. Then the length L of the image of T[y, p] 
under the mapping R satisfies 

L < 2tts ■ max{|i?(C)| : C lies on T[y,p]}. 

Proof. Immediate from Lemma 2.2.1. □ 



2.3 Notes and remarks 

The partial fraction decomposition given in Section 2.1.2, as well as elementary properties of holo- 
morphic functions, can be found in most of the numerous handbooks on complex analysis, e.g. 
Ahlfors (1966), Henrici (1974), Rudin (1974). 

For Stirling's formula, given in Section 2.1.2, see e.g. Henrici (1974); for Cauchy's derivation 
of formula (2.2.2), see Cauchy (1888). 
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In the special case where R(C) is a rational function of order s, with all of its (possible) poles 
at £ = 7, Theorem 2.2.2 is an easy consequence of Bernstein's inequality (see e.g. Edwards (1967)). 
This inequality tells us that, for any such R(Q, we have 

p- max \R'{()\ < 8- max \R(()\. 

IC-7l=P IC~7l=P 

In view of 

L < 2tt P - max \R'{()\, 

C~7l=P 

one arrives immediately at the inequality given by Theorem 2.2.2. 

For the case of general rational fuctions R((), Theorem 2.2.2 was proved in Spijker (1991). 
The arguments used above in the proof of Lemma 2.2.1 are similar to those in that reference. 

A beautiful extension of Theorem 2.2.2, dealing with the image of a circle under the mapping 
R(C) on the Riemann sphere, was given by Wegert & Trefethen (1994). These authors gave a short 
proof of their extension by using a general formula, for the arc length of curves lying on a sphere, 
due to H. Poincare. 
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3 Linear algebra 



3.1 The Jordan canonical form 



The set of all s x s matrices A = (atjk), with complex entries ctjk, is denoted by C s ' s . With the 
usual addition A + B and multiplication 7 • A, for A,E>£ C SyS and 7 g C, the set C s ' s becomes a 
complex vector space of dimension s 2 . 

The s x s identity matrix and zero matrix are denoted by 

/ and O, 

respectively. Further, for all s x s matrices A we define 



A c 



I. 



For any square matrix A we denote the set of its eigenvalues by cr[A]. This set is called the 
spectrum of A. The spectral radius is defined by 

p(A) = max{|A| : A G a[A}}. 

Note that A is regular, i.e. A -1 exists, if and only if G" a [A]. 

Without proof we state the following fundamental theorem. 

Theorem 3.1.1 (Jordan canonical form) Let A be a given sx s matrix. Then there exist a regular 
T G C s ' s and block-diagonal J = diag (J u J 2 , . . . , J r ) G C s ' s with 

A = TJT~ X . 

Here the blocks Jk are bidiagonal matrices of order Sk, 



Jh 



1 

Afc 



1 











1 

Afc 



with Afc G a[A] for k = 1, 2, . . . , r. 



The Jordan matrix J is unique up to permutations of the Jordan blocks Jk- 



Since the order of Jk is Sk we see that, for A G <r[A], the sum 

A fc =A 

is equal to the so-called algebraic multiplicity of A, i.e. the multiplicity of A as a zero of the 
characteristic polynomial of A. Clearly 

si + s 2 H h s r = s. 

We denote the identity matrix of order Sfc by Ik, and define Ek by 

Jk = Afc/fc + E k . 

Further we define s x s matrices i?fc by 

Pk = T diag (O, . . . , 0, I k , 0, . . . , O) T~\ R k = T diag (0, . . . , O, E k , O, . . . , 0) T~\ 

From the above one easily obtains the following theorem (in which Sj k is the Kronecker index 
already introduced in Subsection 1.2.5.2). 
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Theorem 3.1.2 (Spectral representation theorem) Let A be a given s x s matrix. Then there is 
an integer r, with 1 < r < s, such that A can be represented in the form 



A — ^(AfcPfc + R k 



fc=i 



Here \ k £ a[A], PjP k = 5 jk Pj, P k R k = R k P k = R k for 1 < k < r, 1 < j < r, and P x + P 2 + 
■ ■ ■ + P r = P Further, for 1 < k < r, we have P k / O, and there are integers s k > 1 such that 
(R k ) m = O if and only if m> s k . The integers s k satisfy si + s 2 + • • • + s r = s. 



For £ o~[A] we shall often deal with the matrix 

(CI-A)-\ 

It is called the resolvent of A at (. From the spectral representation of A one easily obtains the 
representation of the resolvent as given in the following theorem. 

Theorem 3.1.3 (Representation for the resolvent) Let A G C s ' s and ( cf[A]. Then 

r 

(Ci - A)- 1 = [(C - AfcJ-'P* + (C - \ k )- 2 R k + ■■■ + ((- \ k )- Sk (Rk) Sk - 

k=l 

where X k , P k , R k , s k are as in the spectral representation theorem. 



We denote the complex conjugate of any complex number 7 by 7*, and recall that the 
Hermitian adjoint A* of A = (aj k ) G C s ' s is the s x s matrix whose entry in the j-th row and A:-th 
column equals a£ ■ (for 1 < j < s, 1 < k < s). The matrix A is said to be normal if AA* = A* A, 
and unitary if A A* = I. 

Without proof we state the following important theorem about normal matrices. 

Theorem 3.1.4 (Jordan canonical form of normal matrices) Let A be a given sx s matrix. Then 
A is normal if and only if a Jordan decomposition A = TJT~ l exists with T unitary and all blocks 
J k of order 1. 



3.2 Norms and e-pseudo- eigenvalues 

With C s we denote the vector space of all column vectors 

6 

x= . , with complex £1, £ 2 > • • • , Cs- 
We recall that the Hermitian adjoint x* of x is the row vector given by 

x* = (&,&■■■,£)■ 

Arbitrary norms on C s will be denoted by | ■ |, and the so-called Holder norms by 
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(• 1 1/P 

w P = jEi^-r| for i ^p < °°' 

\x\ p = max \xj\ for p = oo. 

l<j<s 

We thus have |a?| 2 = Vx*x, also called the Euclidean norm of x. 

By a linear functional on C s we mean a linear mapping from C s to C. Such a mapping can 
be represented by a row vector with s complex components. Without proof we state the following 
lemma, which is useful in various practical situations. 

Lemma 3.2.1 (On linear junctionals) Let | • | be an arbitrary norm on C s , and y G C s . Then 
there exists a linear functional F on C s such that 

F(y) = \y\ and \F(x)\ < \x\ for all x in C s . 



Let I • I be an arbitrary vector norm on C s . For any s x s matrix A we write 

P|| = sup \Ax\/\x\, 

the supremum being over all x G C s , The function || • || is easily seen to be a norm on C s ' s , 

the so-called matrix norm induced by the vector norm | • |. We denote the matrix norms, induced 
by the Holder norms, by || • || p . We have 

s 

H^Hoo = max Vlajfcl, ||A||i = ||-4*||oo, H^lh = \/p{A*A). 

1< j < s z — ' 

- J - k=l 

The norm ||^4||2 is also called the spectral norm of A. Clearly, for unitary A we have ||^4||2 = 1- 

In the rest of this Section 3.2, || • || will stand for a matrix norm induced by an arbitrary 
vector norm | • | on C s . 

Since || • || is a norm on the vector space C s ' s , we can measure the 'distance' between two 
s x s matrices A and B by the quantity ||^4 — Accordingly, for given A^ G C s ' s , we shall 
write lim A^ = B and ^2^ =0 A^ = B to denote the situation where B G C s ' s is such that 

k — >oo 

lim \\B — Ak\\ = or lim \\B — ^22=0 ^-k\\ = 0, respectively. We note that the limit concept 

thus defined in C s ' s is independent of the underlying norm || • ||. This is a consequence of the fact 
that for any two norms, say || • || and || • ||', there is a constant a > such that \\A\\' < a ■ \\A\\ (for 
all A G C s ' s ). 

It is easily seen that, for any A, B G C s ' s , we have 

\\AB\\ < \\A\\ \\B\\, and ||/|| = 1. 

From Theorem 3.1.2 (or from the arguments to be given in Example 3.3.7) one can obtain the 
spectral radius formula 

lim p™!! 1 /" = p(A). 

n— >oo 

In the following definition we denote by e an arbitrary, nonnegative real constant. 
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Definition 3.2.2 A is an e-pseudo-eigenvalue of A £ C s ' s if 

(i) there is an £ £ C s ' s with \\E\\ < e such that A G a[A + E]. 

The set of all e-pseudo-eigenvalues of A is called the e-pseudospectrum, and is denoted by cr e L4]. 

□ 

Note that in general the answer to the question of whether A is an e-pseudo-eigenvalue of A 
not ony depends on e and A, but also on the vector norm in C s by which the matrix norm || • || is 
induced. Clearly, for e = 0, the e-pseudospectrum of A is equal to the spectrum a[A]. 

The concept of an e-pseudo-eigenvalue can be related to the following properties. 

(ii) There is an s x s matrix E with = e such that A G cr[A + E]. 
(hi) There is an s x s matrix U with \\U\\ = 1 and \\(A — XI)U\\ < e. 

(iv) There is a vector tie C s with \u\ = 1 and \(A — XI)u\ < e. 

(v) A — XI is singular, or A — XI is regular with \\(A — A/) _1 || _1 < e. 

Theorem 3.2.3 (Characterizations of e-pseudo-eigenvalues) For given e > 0, s > 2, A G C s ' s the 
properties (i)-(v) are equivalent to each other. 



3.3 The Dunford- Taylor integral 

Let tpjk(,C) be holomorphic complex functions on an open subset £1 of the complex plane, where 
1 < j < s, 1 < k < s. We define the mapping from Q into C s ' s , by 

$(0 = (<Pjk(0) for (en. 

For any rectifiable curve T in we define the integral of along T by 
Let || • || be any norm on C s ' s . Then the inequality 



< j r \MO\\ MCI 



holds. It is a consequence of the fact that the norm of a sum does not exceed the sum of the norms 
of the individual terms, and that the integral of $(£) along T can be seen to be a limit of Riemann 
sums. 

Let A be a given s x s matrix, and /(C) holomorphic on 17. Under the assumption that 
cr[A] C we define f(A) by the so-called Dunford- Taylor integral. 



f(A) = ± [ Jj(CW-A)- 1 dC. 



Here T consists of a finite number of rectifiable simple closed curves with positive orientation. 
The interiors fifc of are assumed to be mutually disjoint and to be such that 

a[A] c [jn k , [j (n k ur k ) c a 
k k 

We note that f(A) does not depend on T as long as T satisfies these conditions. This follows by 
applying Cauchy's integral theorem to the entries (pj k (() of <&(£) = /(C) (C^ — A) -1 and using the 
above definition of the integral along T. 

Without proof we state the following basic theorem about the Dunford- Taylor integral. 
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Theorem 3.3.1 Let A G C s ' s and f(Q, ff(C) holomorphic on an open set Q with cr[A] C SI. Then 
(a) /(C • J) = /(C) • I for all ( G fi, 

(bj (a-f + P-g)(A) = a-f(A) + 0-g(A) for all a, (3 G C, 

(cj (/^)(A) = /(A)^(A), 

(dj a[/(A)] = /(<rL4]). 

Assume, additionally: B G C s ' s with Al? = BA; is connected and, for each A G cr[A], the closed 
disk with center A and radius p{B) belongs to ft. Then cr[A + B] C fl and 

oc 1 

(e) f{A + B)= £ m [A)B K 

k=o «! 

Note that the addition and multiplications in the left-hand members of (b), (c) involve scalar 
functions, whereas those in the right-hand members s x s matrices. The right-hand member of (d) 
stands for {/(A) : A G <r[A]}. 

Example 3.3.2 Let R(Q = P(Q/Q(Q be a rational function and A an s x s matrix. We shall 
say that R(A) exists if Q(A) / for all A G cr[A\. In this case the matrix Q(A) is regular by (d), 
and from (c) it can be deduced that 

R(A) = PWiQiA)]- 1 = [0(A)]- 1 P(A). 

Using (b), (c) it can also be shown that the partial fraction decomposition of R(Q (see Theo- 
rem 2.1.1) applies to R(A). If R{A) exists, we have 

R(A) = R (A) + R 1 (A) + --- + R q (A), 

m — n kj 

r (a) = <*i,oA l , R M) = 52<*-iA A ~ Q 1 )' 1 ( for 1 < 3 < q), 

1=0 1=1 

where the coefficients a^j and the integers q, kj are as in Theorem 2.1.1. □ 
Example 3.3.3 Let /(C) = e c . Writing e A = exp (A) = f(A), it follows from (e) that 



OO 



kV 

k=0 



Further, for real t and positively oriented circle T = T[0, p] with p sufficiently large, we have 



2m 



Using the last of these two integral representations it can easily be seen that, when h — > 0, the 
matrix h~ x [e ( ~ t+h ' >A - e tA ] tends to Ae tA . We thus have 

|.( e ^) = Ae tA . 
dt K ' 
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Example 3.3.4 Let A, B be commuting s x s matrices. From (e) we have 



(A + B) n = J2 ( n k ]A k B n ~ k 
fc=0 ^ ' 



We state the following neat theorem without proof. 

Theorem 3.3.5 Let /(£) be holomorphic on an open set O containing the spectrum of the matrix 
A. Let g(£) be holomorphic on f(Q), and define h(Q = g(f(Q) (for ( £ fl). Then h(A) = g(f(A)). 



By substituting the representation for the resolvent, given in Theorem 3.1.3, in the Dunford- 
Taylor integral one arrives after an easy calculation at the following theorem. 

Theorem 3.3.6 (Spectral representation for f(A)) Let A be a given s x s matrix, and /(£) 
holomorphic on an open set Q containing cr[A]. Then 



(3.3.1) f(A) = [/(Afc)Pfc + f'(X k )Rk + j/ 2) (\k)(Rkf 
k=i 



+ ■■■ + ■ 



f^-'HXkXRk) 



Sfc-1 



(s k -l)\- 

Here \k, Pk, Rk, Sk are as in Theorem 3.1.2 (Spectral representation theorem). 

Using the notations of Theorem 3.1.1 (Jordan canonical form), we can rewrite formula (3.3.1) 



as 



(3.3.2a) 

where 

(3.3.26) 



f(A) = TDT~\ 



D = diag(£>i,D 2 ,...,AO 



is a block-diagonal matrix composed of matrices of order Sk- Using the definition of given 
in Section 3.1, we see that allows the elegant representation 



(3.3.3c) 



3=0 



^From this representation it follows that is an upper triangular matrix of the form 



(3.3.3d) 



/(A*) 



f {1 \Xk) / (2) (A.) f^-'HXk) 



1! 



f(X k ) 



2! 

/ (1) (A fc ) 
1! 



(s k ~ 1)! 



7 (2) (A fc ) 
2! 

/ (1) (A fc ) 
1! 

/(A fc ) 
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Example 3.3.7. We illustrate Theorem 3.3.6 by using it in a short proof of the spectral radius 
formula given in Section 3.2. For any integer n > 1 and /(£) = C, 71 -, we have by Theorem 3.3.1(c) 

A n = f(A), 

and by Theorem 3.3.6 

In view of the last two equalities there is a constant 7 such that for n = 1, 2, 3, . . . 

\\A n \\ < jn p [p(A)] n , with p = max (s fe - 1). 

k 

Since ||A n || > p(A n ) = [p{A)] n there follows 

p(A) < p™!! 1 /" < (^n p ) 1/n p(A), 
from which we obtain the spectral radius formula by letting n — > 00. 



3.4 The logarithmic norm 

Throughout this section | • | denotes an arbitrary norm on C s , and || • || stands for the corresponding 
induced matrix norm. 

Consider the initial value problem 

U'{t) = AU(t), U(0) = v, 

where v is a given vector in C s with norm \v\ = 1. ^From the relation ^(e tj4 ) = Ae tA (obtained 
in Example 3.3.3), the solution U(t) G C s can be seen to be equal to U(t) = e tA v (for t G R). The 
logarithmic norm of A, to be defined below, is useful in estimating the norm of e tA v. 
For real t with t 7^ and v, w G C s we define the difference quotient 

m(u, to; t) = t — 1 + tw \ — \v\]. 

The following simple lemma is quite useful. It can be proved by elementary considerations. 

Lemma 3.4.1 For fixed vectors v, w the quantity m(v, w; t) is isotone with respect to t € R\ {0}. 
Further — \ w \ < m(v,w;t) < \w\. 

In view of this lemma we can define the (finite) quantities 

m_(u, w) = \\m. m(v,w;t), m + (v,w) = ]im^m(v,w;t). 

For A e C s ' s we introduce the notations 

m_(A) = sup m-(v,Av), m + (A) = sup m + (v, Av). 

\v\=l \v\ = l 

Further, for A G C s ' s and t 7^ 0, we define 

p{A;t) = sup m(v, Av;t). 

\v\=l 

Using the definition of || • ||, in terms of the norm | • | on C s , it can be seen that 
p(A;t) = t- 1 \\\(I + tA)- 1 \\~ 1 -l\ fort<0, \t\p(A) < 1, 

and 



p{A-t) = t~ 1 



\l + tA\\ - 1 



for t > 0. 



Theorem 3.4.2 (Representations of the logarithmic norm) We have 

\im_ p.(A;t) = m-(A) = m + {A) = \mi + p{A-t). 
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Any of the four quantities in the theorem is called the logarithmic norm of A, and will be 
denoted by 

V(A). 

Theorem 3.4.3 (Characterization of the logarithmic norm) The logarithmic norm of A is equal 
to the smallest constant 00 such that 

\\e tA \\ <e tuj for all t > 0. 

Proof. 1. Prom elementary numerical analysis it is known that the explicit Euler method, for the 
numerical solution of the initial value problem U'(t) = AU(t), U(0) = v, is convergent. Hence 
(J + ±A) n v tends to e tA v as n — > 00. Therefore (I + —A) n tends to e tA (for n — > 00). We have, for 
t > 0, 

||( 7 + l A) n\\ < {1 + 1 KA . i)}" = {1 + l [fi{A) + en]} « > 

with e n — ► for n — ► 00. Consequently ||e* A || < e f ^ A ^ for t > 0. 
2. Suppose w is such that ||e* A || < e* w (for all t > 0). 

Let v £ C s with = 1, and let U(t) denote the solution to the initial value problem 

U'(t) = AU(t), 17(0) = v. 

For t > we have 

m(v, Av; t) = t' 1 [\v + tAv\ - l] = i _1 [|E/(t)| + e(t) - l] , 

where 

lim t _1 • e(t) = 0. 

Hence 

m + (v,Av) = lim mfv.Aujt) < lim [t _1 (e* w - 1) + t _1 e(t)l = w. 
t-»o+ t->o+ L J 

Since /u(A) = m+(A), it follows that <u. □ 

Theorem 3.4.4 (Properties of the logarithmic norm) Let A,B<E C s ' s and 7 G C. Then the 
following relations are valid. 

(i) -\\A\\ < n(A) < \\A\\, 

(ii) n(I) = l, n(O) = 0, l x(A + 1 -I) =fx(A) + Re 7, 

(iii) ^(7 ■ A) = 7 ■ //(-A) provided 7 > (positive homogeneity), 

(iv) + 5) < jLt(^) + (sub-additivity), 

(v) - yii(-B) I < ||A - B|| (continuityj, 

(vi) max{Re A: AG a[A]} < fi(A). 



We denote the logarithmic norm of A corresponding to the norm | • \ p on C s by 

H P (A). 

Theorem 3.4.5 (Expressions for n P (A)) For arbitrary A G C s ' s we have 

Hoo{A) = max (Re (a j:j ) + ^ \a jk \), l^i(A) =/x QO (A*), 

112(A) = max Re (v* Av) = max{ A : A G cr[\(A + A*)]}. 

For normal matrices A we have 

112(A) = max{ Re A : AG op]}. 
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Definition 3.4.6 The matrix A £ C s,s is said to satisfy a circle condition, with respect to the disk 

\\A - 7 • I\\ < p. 



In some of the following it will be important to check whether A satisfies a circle condition, 
and to interprete such a condition in terms of \\A\\ and n(A). The following theorem is useful in 
this context. 

Theorem 3.4.7 {On circle conditions) Let A G C s ' s and a > to. 

a) If\\A+\{a-u)I\\<\{a + u)then\\A\\<a, p{A)<uj. 

b) Assume all diagonal elements a 33 of A = (a 3 k) are real, and p = 1 or p = oo. Then 

\\A + — w)7||p < |(a + u) if and only if \\A\\ p < a, fJ> p (A) < uj. 

Proof, a) Since 

P|| - \{a - uj) < \\A + \{a - uj)I\\ <\{a + uj), 

we have \\A\\ < a. 

Writing 7 = i(a — uj) we obtain 

H(A) = n{A + 7/) - 7 < \\A + 7/H - 7 < (7 + uj) - 7, 

and therefore < uj. 

b) Since for any matrix A we have ||A||i = H^*^, Hi{A) = ^^(A*), it is sufficient to 
prove the statement for p = 00. 

Let || ^4 1| 00 < ck, Moo(^4) < For 1 < j < s we introduce 

°i = J2\ a Jk\- 

From the expressions, given in the above for ||^4||oo, Moo {A), we see that 

\a 33 \ + a 3 <a, • Oj < 

In order to prove ||(A + \{a — w)/||oo < \{ol + w) we only have to show that 

for 1 < j < s. 

Fix j, and define 6 by 

= if a^- + \{a - uj) > 0, 
6 = 1 if + \{a - uj) < 0. 

We have 

\ a jj + U a - UJ )\ + ° j = 6{-a 33 -\{a-uj) + o-j} + (1 - 6){a 33 + \{a-uj) + a 3 } 

< 6{\a 33 \ - \{a -uj) + a 3 } + (1 - 6){a 33 + \{a - uj) + a 3 } 

< 9{a _ 1 (Q _ w)} + (1 _ + 1 (Q _ w)} 

= |(a + w). 

□ 
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3.5 The M-numerical range 

The classical numerical range of a given matrix A £ C s,s is denned as the set of complex numbers 

{x*Ax : iGC s with x*x = 1}. 

It is known to be a convex set containing all eigenvalues of A. Below we consider a generalization 
which will be useful in the following sections. 

Let | • | be an arbitrary norm on C s , and M > 1. For a given matrix A G C s ' s we define the 
disk D[y,p] to be suitable if 

(3.5.1) \\(A-~/I) k \\ < Mp k for A; = 1,2,3, 

Using the spectral radius formula (Section 3.2) we see that, for any suitable disk D[y,p\, the 
inclusion a[A — 7/] C D[0, p] is valid. This implies that 

(3.5.2) a[A}cD[ 7 ,p}. 

In view of (3.5.2), the best enclosure of cr[A] obtainable from relations of the form (3.5.1) 
equals the intersection of all suitable disks. This is a motivation for the following definition. 

Definition 3.5.1 The M-numerical range of A with respect to the norm | ■ | is the set t[A] given 
by 

r[A] = f]D[ 7 ,p], 

where the intersection is over all disks that are suitable. If we want to express the dependence of 
t[A] on M, we write 

t[A]=t[A,M], 

and the numerical ranges corresponding to the Holder norms | • | p are denoted by 

t p [A,M]. □ 

Theorem 3.5.2 (Basic properties of the M-numerical range) For arbitrary norm \ • \ and constant 
M > 1 the following holds: 

(i) t[A, M] is a closed, bounded, convex subset of the complex plane, 

(ii) t[{ I + OA M] = Co + Ci • t[A, M] for all Co, Ci € C, 

(iii) t[A, M 2 ] C t[A, Mi] for 1 < M 1 < M 2 , 

(iv) conv cr[A] C t[A,M], 

(v) f| t[A,M] = conv a[A}. 

M>1 

Proof. The first four properties follow easily from the above definition. In order to prove property 
(v), choose any 70 G C, po > such that conv a[A] lies in the interior of the disk _D[7 ,po]- Since 
p(A — 70-Z") < po the spectral radius formula (Section 3.2) tells us that 

1 

lim \\(A - 7o^) fc || fc < Po- 

Therefore, there is an M with 

||04 - 7 o/) fc || <M (p ) fc for fc = 1,2,3, ... , 

so that t[A,M ] C £>[7o,Po]- 

Since conv a[A] is equal to the intersection of all disks D[j ,p ] of the above type, we see 
that conv a[A\ must be equal to the intersection of all sets t[A, M] with M > 1, i.e. (v). □ 
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Let W be an arbitrary convex subset of C, and let ( £ C. The distance from ( to W is 
defined by 

d((,W)=mi {|C-el : i^W}. 
If £ belongs to the boundary dW of II 7 and 

Re(e" ie (C - 0) < for a11 C G W, 

where 6 is a real constant, then 6 is said to be a normal direction to W at £. 

We shall deal with the following four conditions on A G C s ' s with respect to W C C. 



(I) r[A,M]cW. 

(II) C-f - A is regular and ||(C/ - A)~ k \\ < M • [d((,W)]~ k for all ( G C \ W and 
fc = 1,2,3,.... 

(III) ||exp [te- w (A - £J)]|| < M for alU > 0, £ G cW, and normal directions (9 to W at f. 

(IV) There is a vector norm | • 1 on C s such that the corresponding 1-numerical range 
satisfies t [A, 1} C W and \x\ < \x\ < M ■ \x\ (for all x G C s ). 



Theorem 3.5.3 (Main theorem on the M-numerical range) Let \ ■ \ be an arbitrary norm on C s , 
and A G C s,s . Let M > 1, and W an arbitrary nonempty, closed, and convex subset of C. Let 
t[A,M] be the M-numerical range of A with respect to \ ■ \. Then the conditions (I)-(IV) are 
equivalent to each other. 

Clearly, t[A, M] is the smallest nonempty, closed and convex set W C C with property (I). 
Therefore, the above theorem reveals three new characterizations of the M-numerical range. We 
see that r[A, M] equals the smallest nonempty, closed and convex set W C C with property (II), 
and the same holds with regard to the properties (III) and (IV). 

We note that, for M = 1, further characterizations of the set t[A, M] are possible. One 
of these is easily obtained by applying the characterization of the logarithmic norm (see Theorem 
3.4.3 with u = 0) to the norm-inequality in (III). We see that, for M = 1 and | • |, A, W as in the 
above theorem, the conditions (I)-(IV) become equivalent to the condition 

(V) /i(e _i M) < Re(e" ie £) for all £ G dW, and normal directions 6 to W at £. 

Theorem 3.5.4 (Expressions for t p [A,1\) For arbitrary A = (ajk) G C s ' s let A' = (a' jk ) with 
a'- k = otkj, and let 

Dj = D\ujj,Gj\, a j = ^2\a jk \ for j = 1, 2, . . . , s. 

We have 

s 

^ [A, 1] = conv (J Dj , n [A, 1] = r OT [A' , 1] , 

i=i 

t 2 [A, 1] = {x*Ax : x G C s with x*x = 1}. 

For normal A one has 

t 2 [A, 1] = conv o[A\. 
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Using the above expressions for t p [A, 1] with p = l,oo one can establish a simple relation 
between circle conditions and the set t p [A, 1]. We have the following theorem. 

Theorem 3.5.5 (Relation between circle conditions and T p [A, 1] for p = l,oo) Let D[y, p] be an 
arbitrary disk, A = {oijk) € C s ' s and p = 1 or oo. Then A satisfies a circle condition with respect 
to D[j, p] and \\ • \\ p , if and only if t p [A, 1] C D\y, p\. 

Proof. Since = 11^4'Hoo, T\[A, 1] = TqoLA',1], where A' is as in Theorem 3.5.4, we only have 

to consider the case p = oo. 

If A satisfies a circle condition with respect to D[y,p\, we have Too [.A, 1] C D[y,p] by the 
definition of the 1-numerical range. 

Conversely, assume t^A, 1] C D[y,p\. In view of Theorem 3.5.4, and with the notations of 
that theorem, there follows 

D[aj3,aj] C D[j,p], 

and consequently \ctjj — j\ + <jj < p (for j = 1, 2, . . . , s). Hence \\A — 7-f ||oo < P- □ 
3.6 Linear algebra concepts and semi-discretization 

We shall illustrate some of the above concepts and theorems by applying them to the matrices A 
originating in the process of semi-discretization described in Section 1.2.5. 
Consider problem (1.1.6) and assume 

(3.6.1) a(x) > 0, b(x) < 0, c(x) < 0, and b(x) is antitone . 

We shall deal with the matrix A defined by (1.2.7b), (1.2.8), (1.2.9a), (1.2.10a). Assume, with 
regard to (1.2.9a), that 

(3.6.2) all upwind parameters satisfy £\ = 1. 
In this situation the matrix A = (ctjk) is tridiagonal; we have 

(3.6.3) a j,j+i = ^~ 2a j+\ (f° r 1 — 3 — s ~ 1)' 

ctjj-i = S~ 2 Oj_i — 5~ 1 b j _i (for 2 < j < s — 1), 
= 5 ~ 2 + a s+ i] - <J _1 & s _i (for j = s), 

From (3.6.1) we have, for 2 < j < s, the inequalities 

and 

an + y~]|oiifc| = — 5~ 2 ai + 5~ 1 b3 + Cj < 0. 
fc/i 

Therefore, by Theorem 3.4.5, 



(3.6.4) 



fioo(A) < 0. 
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By virtue of Theorem 3.4.3 there follows 

||exp (tA)Hoo < 1 for all t > 0. 

The last conclusion can be interpreted as a stability result for the semi- discrete problem 
(1.2.7a). Suppose U (t) is a solution to (1.2.7a) with uo replaced by uq. Then 

U(t) - U(t) =exp (tA){u -u ), 

and since the above matrix norm of exp (tA) does not exceed 1, we have 

\U(t)-U(t)\ < \u -u \ 



I oo 



This means that any perturbation in the initial vector uq is not amplified, when the time t increases 
and differences between solutions are measured in the maximum norm | • | oo • This stability result for 
(1.2.7a) is a semi-discrete analogue of a stability property for (1.1.6) which is plausible from (3.6.1) 
when viewing a(x), b(x), c(x) as corresponding to diffusion, convection and reaction, respectively. 

An essential point in the above derivation of (3.6.4) is the property of A that all of its off- 
diagonal entries are nonnegative. If (3.6.2) is not fulfilled, then this property need not be present, 
and accordingly (3.6.4) can be violated. This explains once more why upwind approximations are 
advantageous. 

In the situation (3.6.2) one can estimate ^^{A), similarly as above, also for arbitrary a(x) > 
and b(x), c(x) not necessarily satisfying (3.6.1). In this case (3.6.4) is replaced by an inequality 

(3.6.5) Hoo(A) < lo 

with a constant u>, depending only on the functions b(x), c(x). 

The above estimates for ^^(A) can still be established for values e\ somewhat smaller than 
1. In fact, if 

(3.6.6) e x >l- for A=-,-,...,— , 

then all off-diagonal entries of A are still nonnegative, and the same estimates are valid as for the 
fully upwind approximations defined by (3.6.2). The ratio 

a\ 

is called the (mesh-)Peclet number at the location x\ = A • 5. If diffusion dominates convection, we 
have small Peclet numbers, and (3.6.6) may be fulfilled with £\ = 0. But, for strongly convection 
dominated problems the Peclet number is large, and (3.6.6) forces e\ to be close to 1. 

We turn again to the situation (3.6.1), (3.6.2). From the general expression for ||-A||oo (cf. 
Section 3.2) we see that the tridiagonal matrix A defined by (3.6.3) satisfies 

(3.6.7) PIU < a, 
with 

(3.6.8) a = 4<T 2 |a| 00 + 2S~ 1 \b\ 00 + (cU- 

By applying Theorem 3.4.7 (part b), with u = and p = oo, it follows that A satisfies the circle 
condition 

(3-6.9) IIA+Hju^H. 

This inequality implies that the 1-numerical range ^[^4, 1] is contained in the disk D[—^,^]. More 
refined information about the shape of t^A, 1] can be obtained by using Theorem 3.5.4. 
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3.7 Notes and Remarks 

A proof of Theorem 3.1.1 (Jordan canonical form) and Theorem 3.1.4 (Jordan canonical form of 
normal matrices) can be found e.g. in Horn & Johnson (1990). The representation Theorems 3.1.2, 
3.1.3 can be derived directly, without using Theorem 3.1.1, by using complex integration theory, 
see e.g. Kato (1976). 

For further discussion of norms in C s and C s ' s we refer to Horn &i Johnson (1990). Lemma 

3.2.1 (Linear functionals) can be viewed as a corollary to the so-called Hahn-Banach theorem from 
functional analysis, cf. Rudin (1973). For a direct derivation of Lemma 3.2.1, in the context of the 
finite dimensional space C s , see e.g. Horn & Johnson (1990). 

The spectral radius formula (Section 3.2) is a special case of a general result valid in complex 
Banach algebras, cf. Rudin (1973). For further discussion of the spectral radius formula, for 
matrices A G C s,s , we refer to Horn & Johnson (1990). 

The concept of the e-pseudospectrum was introduced and studied a.o. by Landau (1975), 
Reddy & Trefethen (1990, 1992), Varah (1979). Very interesting properties of the e-pseudospectrum, 
as well as a wealth of applications, were discovered by L.N. Trefethen, see Trefethen (1996). The 
focus in the works just mentioned is on the (weighted) spectral norms \\A\\ 2 for A G C s ' s . The 
e-pseudospectrum for the case of arbitrary norms \\A\\ in C s,s was discussed in Dorsselaer, Kraaije- 
vanger & Spijker (1993). Theorem 3.2.3 (Characterizations of the e-pseudo-eigenvalues) has been 
taken from that paper. 

For additional reading about the Dunford- Taylor integral, see Conway (1985), Dowson 
(1978), Dunford & Schwartz (1958), Kato (1976). A direct and neat proof of the Theorems 3.3.1, 
3.3.5, valid for the general case of linear operators in a Banach space, is given in Dowson(1978). 
We note that the parts (a)-(d) of Theorem 3.3.1 can also be proved easily by using Theorem 3.3.6. 
For further discussion of the latter theorem, see e.g. Dunford & Schwartz (1958), Kato (1976). 

The logarithmic norm was introduced independently by Dahlquist (1959) and Lozinskij 
(1958). Subsequently, various properties of the logarithmic norm were established by Desoer & 
Haneda (1972) and Strom (1975). The last equality in Theorem 3.4.2 (m + (A) = \im fx{A;t)) and 

the Theorems 3.4.3, 3.4.4, 3.4.5 can be found, with complete proofs, in the above four references. 
Lemma 3.4.1 has been taken from Martin (1976, p. 37). A proof of the second equality in Theorem 

3.4.2 (m_ (A) = m + (A)) can also be found in Martin (1976, p. 246) — even for the general situation 

of operators acting on a Banach space. The fact that lim a(A; t) = lim a(A; t) was proved in 

t->o- t^o+ 

Dorsselaer & Spijker (1994). Theorem 3.4.7 (Circle conditions) has been taken from Spijker (1985). 

There exists a vast literature on the classical numerical range and closely related issues, see 
Horn <fe Johnson (1994). For M = 1, the M-numerical range (Definition 3.5.1) can be seen to be 
equal to the so-called algebra numerical range, well known in some parts of functional analysis, see 
Bonsall & Duncan (1973, p. 42) and (1980, p. 35). The M-numerical range, for general M > 1, was 
introduced in Lenferink & Spijker (1990). In that paper Theorem 3.5.2 and part of Theorem 3.5.3 
(equivalence of (I), (II) and (III)) was proved. The equivalence of (I) and (IV) was obtained in 
Spijker (1993). The expressions for t p [A, 1], given in Theorem 3.5.4, have been known for a long 
time, cf. Lenferink & Spijker (1990), Spijker (1993) for corresponding historical remarks. 
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4 The problem of stability in the numerical solution 
of differential equations 

4.1 Linear stability analysis 

We shall deal with step- by-step methods for the numerical solution of linear differential equations. 
Both initial-boundary value problems in partial differential equations and initial value problems in 
ordinary differential equations will be included in our considerations. 

A crucial question in the step-by-step solution of such problems is whether the method 
will behave stably or not. Here we use the term stability to designate that any numerical errors, 
introduced at some stage of the calculations, are propagated in a mild fashion - i.e. do not blow 
up in the subsequent steps of the method. 

Classical tools to assess the stability a priori, in the numerical solution of partial differential 
equations, include Fourier transformation and the corresponding famous Von Neumann condition 
for stability. Further tools of recognized merit for assessing stability, in the solution of ordinary 
differential equations, comprise so-called stability regions in the complex plane. Since the mid sixties 
these stability regions have been studied extensively; numerous papers have appeared dealing with 
the shape and various peculiarities of these regions. 

However, the above tools are based on the behaviour that the numerical method would have 
when applied to quite simple test problems. Accordingly, in the case of partial differential equations 
Fourier transformation provides a straightforward and reliable stability criterion primarily only for 
certain numerical methods applied to pure initial value problems in linear differential equations with 
constant coefficients. In many cases of practical interest, Fourier transformation is not relevant to 
analysing stability: e.g. for pseudo-spectral methods applied to initial-boundary value problems, 
for finite difference or finite element methods related to highly irregular grids, and for methods 
applied to equations with strongly varying coefficients. Similarly, in the case of ordinary differential 
equations, stability regions are primarily relevant only to numerical methods when applied to scalar 
equations 



with given complex constant A. 

Clearly, rigorous stability criteria with a wider scope than the simple classical test equations 
are important - both from a practical and a theoretical point of view. It is equally important to 
know to what extent stability regions can be relied upon in assessing stability in the numerical 
solution of differential equations more general than (4.1.1). In the following we shall discuss var- 
ious theories which are relevant to these two questions. No essential use will be made of Fourier 
transformation. 

4.2 Stability and power boundedness 

In the following we shall deal with numerical processes of the form 

(4.2.1) u n = Bu n _ 1 + r n for n= 1,2,3,..., 

with a given square matrix B of order s > 1 and given vectors r n £ C s . The s-dimensional vectors 
u n are computed sequentially from (4.2.1) starting from a given vector uq G C s . 




U'(t) = XU(t) for t>0, 
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Processes of the form (4.2.1) occur in the numerical solution of linear initial value problems 
that are essentially more general than the simple classical test problems mentioned above. The 
vectors u n provide numerical approximations to the true solution of the initial(-boundary) value 
problem under consideration. 

As an illustration to (4.2.1) we consider again the one-dimensional diffusion-convection- 
reaction problem (1.1.6) of Section 1.1, 

(4.2.2) ^ t < x ^) = ^M*)^^'*)] + -^[b(xHx,t)} + c{x)u{x,t) + d{x), 

d 

u(0, t) = g(t), q^u(1, t) = 0, u{x, 0) = f(x), where < x < 1, t > 0. 

We choose positive increments At = h, Ax = 5 = 1/s and consider the approximation of u(x,t) 
at x = Xj = j5, t = nh. We denote these approximations by u 7 -. The following finite difference 
scheme has been constructed by applying the 2-stage 6>-Runge-Kutta method (an = a\ 2 = 0, 
a 2 i = 6i = 1 - 9, a 22 = b 2 = 6) to (1.2.7), (1.2.8), (1.2.9a) (with all e x = 1), (1.2.10a). 

h-^u] - u]- 1 ) =5~ 2 {a j+1/2 [6u] +1 + (1 - 9)u]~l] - (a j+1/2 + a 3 _ 1/2 ) [6u]+ 
+ (i - e)^- 1 ] + a,_ 1/2 [Ou^ + (1 - 9)u-zl] } 
+0" 1 {b j+1/2 [0u? + (l - e)^" 1 ] -6,_ 1/2 [^_ 1 + (1 - 0)^] } + 
+ Cj + (1 - fl)^" 1 ] +dj, 



,n-l 



5 ((n-l)n), 



n-1 
l s-l > 



where j = 1,2, ... ,s and n = 1, 2, 3, . . .. In the above 6 denotes a parameter, with < < 1, 
specifying the numerical process, and xa = A<5, oa = o(xa), »a = b(x\), c\ = c(x\), d\ = d(x\). 
Defining vectors u n by 

fv%\ ( u(x 1 ,nh)\ 



u(x 2 ,nh) 



\u(x s ,nh) j 



one easily verifies that the u n satisfy a relation of the form (4.2.1). Here 

(4.2.3) B = {I + {1 - 6)hA){I - 6hA)- 1 , 

where A = (ctjk) is an s x s tridiagonal matrix with its (nonzero) entries given by 



(4.2.4) a h j+1 = 5~ 2 a jH (for 1 < j < s - 1), 



2 c 

ctjj-i = d~ 2 a j _i - 5~ 1 b j _i (for 2 < j < s - 1), 
"m-i = 5 ~ 2 [°s-i + a s + 1] - 5_l6 s-i ( for J = a) 



-o 2 [a 3 -_ i + a j+ i] + o + Cj (for 1 < j < s). 



Clearly this matrix A equals the matrix given by (1.2.7b), (1.2.8), (1.2.9a) (with all e\ = 1), 
(1.2.10a), and it coincides with (3.6.3). 

Suppose the numerical calculations based on the general process (4.2.1) were performed 
using a perturbed starting vector uq, instead of uq. We would then obtain approximations that 
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we denote by u n . For instance uq may stand for a finite-digit representation in a computer of the 
true Uo, and the u n then stand for the numerical approximations obtained in the presence of the 
rounding error vo = uq — uq. 

In the stability analysis of (4.2.1) a crucial question is whether the difference v n = u n — u n 
(for n > 1) can be bounded suitably in terms of the perturbation v = uq — uq. Since 

v n = u n -u n = [Bu n _i + r n ] - [Bu n -i + r n ] 

we have 

v n = Bv n - U 

and consequently 

v n = B n v . 

The last expression makes clear that a central issue in our stability analysis is the question 
of whether given matrices have powers that are uniformly bounded. Accordingly, in the following 
we focus, for an arbitrary s x s matrix B, on the stability property 

(4.2.5) ||5 n ||<M for n = 0,1,2,..., 

where M is a positive constant. Throughout this Chapter 4 we denote by || • || the spectral norm 
(see Section 3.2). 



4.3 Power boundedness and the eigenvalue criterion 

For any given matrix B one can easily deduce a criterion for the existence of an Mq with property 
(4.2.5). In view of the Theorems 3.3.1, 3.3.6 we have 



(A/, )" T\ + ( ") (A/,)" '/,', + (\ k ) n - 2 {R k f + ■■■+( k n _ x ) (\k) n - Sk+1 {R k ) Sk - 1 



k=i 

(cf. Example 3.3.7). Here X k , Pk and R k are related to the Jordan decomposition of our matrix B 
similarly as they were related to the decomposition of A in Section 3.1. 

From this representation for B n we easily conclude that ||.B n || remains bounded when n — > oo 
if and only if 

|Afc| < 1, and s k = 1 when |Afc| = 1. 
We thus arrive at the following theorem. 

Theorem 4.3.1 (The eigenvalue criterion) For a given matrix B there exists a constant M with 
property (4.2.5) if and only if 



(4.3.1) 



all eigenvalues A of B have a modulus |A| < 1, and any Jordan 
block corresponding to an eigenvalue A with |A| = 1 has order 1 . 
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However, in the stability analysis of numerical processes one is often interested in property 
(4.2.5) for all B belonging to some infinite family T of matrices. The crucial question then is of 
whether a single finite Mq exists such that (4.2.5) holds simultaneously for all B belonging to T . 
In this situation, (4.3.1) may only provide a condition that is necessary (and not sufficient) for such 
an M to exist. 

For instance, in the example of Section 4.2 one can only expect great accuracy in the ap- 
proximations «™ to u(xj,nh) when 5 (and h) become very small. Accordingly one is primarily 
interested in bounds for B n that are valid uniformly for all B of the form (4.2.3), (4.2.4) with 
arbitrarily small 6 = 1/s. 

An instructive counterexample, illustrating the fact that the criterion (4.3.1) can be mis- 
leading for the case of families J 7 , is provided by the s x s bidiagonal matrices 



(4.3.2) 



B 



(-1/2 

3/2 -1/2 

: 3/2 

V 







3/2 



\ 






-1/2/ 



Matrices of the form (4.3.2) may be thought of as arising in the numerical solution of the 
initial-boundary value problem 

d d 

—u(x,t) = ——u(x,t), u(0,t) = 0, u(x,0) = fix), where < x < 1, t > 0. 
ot ox 

This problem is the same as (1.1.5) with b(x) = —1, c(x,t) = 0. Application of the (forward) Euler 
method to a fully upwind semi-discrete approximation, based on (1.2.1b), yields the numerical 
process 



hr x (u r - 



u 



n-l\ 



.71-1 



= 0, u° = f(j/s). 



Here At = h > 0, Ax = 5 = 1/s, and n" approximates u(jS,nh) for j = 1,2, ... ,s and n = 
1,2,3,.... Clearly, with the choice hs = 3/2 this numerical process can be written in the form 
(4.2.1) with a matrix B as in (4.3.2). 

For each s > 1 the matrix (4.3.2) satisfies the eigenvalue condition (4.3.1). 

Defining the s x s shift matrix E by 



(4.3.3) 



E 



1 


Vo 



1 o/ 



we have from (4.3.2) the expression 
Consequently (cf. Example 3.3.4), 



b—\i + \b. 



Bn = J2 Q(-V2) n ~ fc (3/2) fc £ fc . 
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Defining x to be the s-dimensional vector whose j-th component equals £j = (— 1) J , and denoting 
the j-th component of y = B n x by r\j we easily obtain, from the above expression for B n , 

\7 lj \ = J2('^\(l/2) n - k (3/2) k =2 n for n + l<j<s. 

For s > n we thus have 




Since 




the spectral norm of B n satisfies ||-B n || > y/l — n/s • 2™. Denoting the s x s matrix B by B s we 
thus have 

\\(B 2n ) n \\>2 n - 1 / 2 for n = 1,2,3,... 

Clearly, no M can exist such that (4.2.5) is valid for all B belonging to T = {B s : s = 1, 2, 3, . . .}. 

It should be noted that in some special cases the eigenvalue criterion can be reliable. For 
normal matrices B we have, with the notation of Theorem 3.1.4 applied to the matrix B, 

B n = (TJT -1 ) n = T diag(A?, A£, . . . , A")T _1 , 

so that (4.3.1) implies ||5 n || < ||T|| • ||T _1 ||. Since T~ l = T* we have ||T|| = \\T~ l \\ = 1, and we 
arrive at the following theorem. 

Theorem 4.3.2 {The eigenvalue criterion for normal matrices) Let the matrix B be normal, and 
|| • || denote the spectral norm. Then the following three statements are equivalent to each other. 

(i) There is an M such that (4.2.5) is fulfilled. 

(ii) \\B n \\ < 1 forn = 1,2,3,.... 

(in) All eigenvalues A of B have a modulus |A| < 1. 

But, in general the matrices B in (4.2.1) are not normal, and one has to look for conditions 
different from (4.3.1). In the following chapters we shall deal with conditions which are still reliable 
for the case of non normal matrices. 



4.4 Notes and remarks 

For a discussion of Fourier transformations and the corresponding Von Neumann condition for 
analysing stability we refer to the classical work Richtmyer & Morton (1967) and to Strikwerda 
(1989). 

Already in the pioneering work by F. John (1952) the scope of Fourier transformation was 
widened in that it was used in deriving sufficient conditions for stability in the numerical solution 
of linear partial differential equations with mildly variable coefficients. For subsequent related 
work, relevant to equations with variable coefficients and also to initial-boundary value problems, 
the reader may consult Gustafsson, Kreiss & Sundstrom (1972), Kreiss (1966), Meis & Marcowitz 
(1981), Richtmyer & Morton (1967), Strikwerda (1989), Thomee (1990) and the references therein. 
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Stability regions, related to numerical methods for ordinary differential equations, are dis- 
cussed extensively in the excellent works by Butcher (1987) and Hairer & Wanner (1996). 

The fact that the eigenvalue criterion (4.3.1) can be a misleading guide to stability was 
already known in the sixties, see e.g. Parter (1962). A related, but stronger, necessary requirement 
for stability is the so-called Godunov-Ryabenkii stability condition, a discussion of which can be 
found e.g. in Morton (1980), Richtmyer k Morton (1967), Thomee (1990). The latter condition is 
not satisfied in example (4.3.2). 

The counterexample (4.3.2) is similar to examples in Kreiss (1990), Reddy & Trefethen 
(1992), Richtmyer & Morton (1967), Spijker (1985). Further examples of instability under the 
eigenvalue condition (4.3.1) can be found in Griffiths, Christie & Mitchell (1980), Kraaijevanger, 
Lenferink & Spijker (1987), Lenferink Sz Spijker (1991b). The matrices B in these references have 
s different eigenvalues A with |A| < 1, and occur in the numerical solution of problems of the form 
(4.2.2). See Reddy & Trefethen (1992), Trefethen (1988) for related counterexamples in spectral 
methods. 
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5 Stability estimates under resolvent conditions 
on the numerical solution operator B 

5.1 Power boundedness and the Kreiss resolvent condition 

By || • || we denote, in this Section 5.1, the spectral norm. 

In the early sixties H. -O. Kreiss established an important theorem, nowadays called the 
Kreiss matrix theorem, dealing with the uniform boundedness of ||-B n ||. Here n = 1, 2, 3, . . . and B 
belongs to a given, possibly infinite, family T of s x s matrices (with fixed s > 1). The theorem 
gives a series of conditions which are equivalent to the existence of a constant M such that 

(5.1.1) ||5 n ||<M for n = 0,1,2, ... , 

simultaneously for all B G T ' . 

One of the conditions in the Kreiss theorem involves the resolvent ((I — B)~ l of B, and 
amounts to the requirement that 

(5.1.2) (I- B is invertible, and ||(C/ - £) _1 || < M 1 • (|C| - 1)~\ 
for all complex numbers £ G C \ D. 

Here Mi is a positive constant, and D denotes the closed unit disk {( : ( G C with \(\ < 1}. We 
shall refer to (5.1.2) as the Kreiss resolvent condition. In many cases of practical interest it is easier 
to verify (5.1.2) than (5.1.1). 

Let B be a given s x s matrix satisfying (5.1.1). Then, by Theorem 4.3.1 (Eigenvalue 
criterion), all eigenvalues of B lie in D. For ( € C\D, the matrix (I — B is thus invertible and, in 
view of Theorem 3.3.1, 

oo 
n=0 

It follows that 

oo oo 

m-B)-'\\ < ^ici—mi^ii < mo^ki— 1 = Mo-aci-ir 1 . 

n=0 n=0 

We see that (5.1.1) implies (5.1.2) with M x = M . 

The Kreiss matrix theorem asserts that, conversely, (5.1.2) implies (5.1.1) with M depending 
only on Mi and on the dimension s, but otherwise independent of the matrix B. 

The Kreiss theorem has often been used in the stability analysis of numerical methods for 
solving initial value problems for partial differential equations. In the classical situation, considered 
in the sixties, the matrices B are obtained by Fourier transformation of the numerical solution 
operators, and they stand essentially for the so-called amplification matrices. These matrices are 
of a fixed finite order s. On the other hand, the implication of (5.1.1) by (5.1.2) can also be used 
without Fourier transformation, with B standing for the numerical solution operator in (4.2.1). In 
this situation we may be dealing with a family of matrices B of finite - but not uniformly bounded 
- order s. Therefore, of particular interest is the dependence of the stability constant M in (5.1.1) 
on the dimension s. 

Since the work of Kreiss many authors studied the size of (the optimal) M as a function of 
Mi and s, and eventually in the ninetees some open problems in this field were settled. Moreover, 
the implication of (5.1.1) by (5.1.2) as discussed above was generalized in several directions. More 
general norms than the spectral norm were dealt with and the resolvent condition (5.1.2) was 
adapted to domains different from the unit disk D. In the following we shall discuss some of the 
results just mentioned as well as closely related ones. 
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5.2 Stability estimates for arbitrary Mi > 1 and arbitrary norms 

In this section we consider the relation between (5.1.1) and (5.1.2) for the case where || • || is a 
matrix norm on C s ' s induced by an arbitrary vector norm | • | on C s . Throughout this section || • || 
denotes such a matrix norm. 

Theorem 5.2.1 Let s > 1, B G C s > s . 

(a) If (5.1.1) holds for some M , then (5.1.2) holds with M x = M ; 

(b) If (5.1.2) holds for some Mi, then 

(5.2.1) < (l + l/n) n min{n + l,s}Mi < e • min{n + 1, s} M x for n = 1,2,3,.... 

Proof. 1. The proof of (a) is the same as the proof in Section 5.1 for the spectral norm. 

In order to prove (b) we consider an arbitrary but fixed n > 1 and B satisfying (5.1.2). In 
view of Theorem 3.3.1(c) we can express the n-th power B n as a Dunford- Taylor integral 

where T is any positively oriented circle |£| = 1 + e with e > 0. Taking norms and applying (5.1.2) 
we obtain 

||£ n || < ^- f(l + e) n Mie- 1 \d(\ = (l + e- 1 ){l + e) n Mi. 
By choosing e = 1/n there follows 

||# n || < (l + l/n) n (n + l)Mi. 
Clearly, the proof of (b) is complete if we can show 

(5.2.2) \\B n \\ < (1 + l/n) n sMi for n= 1,2,3,... . 

2. We can regard C s ' s as the space C*, with t = s 2 , and the norm || • || as a vector norm on 
C*. Applying Lemma 3.2.1 (Lineair functionals) , with y = B n and s replaced by t, we see that a 
linear F : C s ' s — ► C exists with 

(5.2.3) |F(A)| < P|| for all sxs matrices A, 

(5.2.4) F(B n ) = \\B n \\. 

A combination of (5.2.4) and the above integral representation for B n yields 

P n ll = ^r^C n ^(CK, 

where R(Q = F(((I — B) -1 ). Integrating by parts we obtain, with e = 1/n, 

< 5 - 2 - 5 > ii B "» = 2^Ti)//'" +lR '(^ £ ^(i + i/nr/woiuci. 

3. Let Ejk stand for the sxs matrix with entry in the j-th row and k-th column equal to 
1, and all other entries 0. Denoting the entries of the matrix ((I — B) -1 by rjk(C) we thus have 

(CJ-B)-^ J> ifc (C)i? ifc , 



39 



and therefore also 

R(() = ^2r jk (C)F(E jk ). 

j,k 

We recall that a rational function is of order s if its numerator and denominator are polynomials of 
a degree not exceeding s. By Cramer's rule, the rjk(Q are rational functions of order s, and they 
have the same denominator. Hence R(Q is also a rational function of order s. 
By Theorem 2.2.2 we have 

(5.2.6) f \R'(()\ |dC|<27rsmax|B(C)|. 

Jr r 

The proof of (5.2.2) now easily follows by a combination of (5.2.5), (5.2.6), (5.2.3) and (5.1.2). □ 

In the following theorem we focus on the sharpness of the bound (5.2.1) in the case n = s — 1. 
Theorem 5.2.2 Let s > 2. Then 



1 



(5.2.7) suplWB^W /M^B) : B G C s ' s , M 1 (B) < oo} = (l + 7— y)* s 



where M 1 (B) denotes the smallest Mi such that (5.1.2) holds (we define M\{B) = oo if (5.1.2) is 
not fulfilled for any Mi). 

Proof. Define B G C s ' s by B = jE, where 7 > is large and the s x s matrix E is defined by 
(4.3.3). We have 



M 1 (B)= su V (\Q-l)\\{QI-B)- 1 \\= sup%i 
ICI>i Kl>i Kl 



s-1 



.7=0 



s-l 



<J ft f|ff|=Mf I |B , " 1 ||(l + 0(7- 1 )), 



where 



so that 



Mi = supdci-ijicr'- 1 

ICI>1 



max (1 — x)x J = j J (j + 1) 

0<rE<l 



-i-i 



|| J B s - 1 ||/M 1 ( J B)>l//x s _ 1 +0( 7 - 1 ) 

= (l + j^y) 5 1 s + 0(7" 1 ) (as 7^00). 

It follows that the left-hand member of (5.2.7) is not smaller than the right-hand member. In view 
of (5.2.1) the proof is complete. □ 

Corollary 5.2.3 For each s > 1, let an induced matrix norm \\ ■ \\^ be given on C s ' s . Then there 
exist matrices B s G C s,s for s = 1, 2, 3, . . ., such that Mi(B s ) < 00 and 



(5.2.8) ||(B s ) a - 1 ||W~e5Mi( J B fl ) (as s - 00), 

where M\(B a ) has the same meaning as in Theorem 5.2.2 (with \\ ■ \\ - 



\ {s) )- 
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Proof. Immediate from Theorem 5.2.2. □ 
In view of (5.2.1), the estimate 

(5.2.9) ||S n ||<esMi for n = 0, 1, 2, . . . 

is valid for general induced matrix norms on C s ' s . By virtue of Corollary 5.2.3, this estimate is 
sharp in the sense of (5.2.8). However, it should be emphasized that this does not resolve the 
sharpness question for given fixed Mi, since M 1 (B S ) in (5.2.8) may depend on s. In the next two 
sections we will focus on the situation where Mi is a given fixed number. 



5.3 Improved stability estimates for Mi = 1 
5.3.1 The case of arbitrary norms 

In the special situation where the resolvent condition (5.1.2) holds with Mi = 1, the upper bound 
(5.2.1) can be improved in various ways. First we concentrate on arbitrary induced norms on C s,s . 

Theorem 5.3.1 Let s > 1, B € C s ' s and \\ • \\ denote an arbitrary induced matrix norm on C s,s . 
If (5.1.2) holds with M x = 1, then 



(5.3.1) ||5 n || < n!(e/n) n < y/2ir(n + 1) for n= 1,2,... . 

Proof. Property (5.1.2) with Mi = 1 implies that condition (II) (of Section 3.5) is fulfilled by A = B 
with M = 1, W = {( : \(\ < 1}. By virtue of Theorem 3.5.3 this A must satisfy the corresponding 
condition (III) as well, i.e. 

||exp [te~ w (B - e w I)}\\ < M = 1 for all t > and real 9. 

Consequently, 

||exp (C-B)H < exp (|£|) for all complex (. 

From the power series expansion for exp(A) (Example 3.3.3 in Section 3.3) with A = (B, it 
can be seen that 

B n = ^ [ ^C n - 1 ew (CB)dC 

where T is the positively oriented circle with radius n and center 0. Therefore ||-B n || < n! n~ n e n . 
In view of Stirling's formula (Theorem 2.1.2) we have 

n\ < (n/e) n v / 2^exp [(12n) -1 ], 

from which it can be deduced that n\ (e/n) n < y/2ir(n + l). □ 

As the next theorem shows, the upperbound for ||.B n || in (5.3.1) is not unnecessarily pes- 
simistic. 

Theorem 5.3.2 Let s > 2, and let B = (flij) denote the s x s matrix with 

Pij = 1 (for j = i- 1), and 0^ = (for j^i - 1). 

Then there is a vector norm in C s such that, for the corresponding induced matrix norm \\ ■ \\ on 
C s,s , we have 

(i) IKCI-i?)" 1 !! < (|C| -I)" 1 for all ( € C with |C| > 1; 

(ii) \\B n \\ = n!(e/n) n > V2vr • n forn = 1, 2, . . . , s - 1. 
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Proof. 1. For any column vector x in C s , we denote its components by xq,Xi, . . . ,x s -i, and we 
introduce the polynomial 

x(t) = xq + x\t H h x s _it s_1 . 



We define 



= inf i ^ |a fc |e |A,t : m > 1, a k G C, A fe G C with ^a k e XkB = x(B) I. 
^fc=i fc=i ^ 



Note that the relation 

m 

]>> fc e A * B =x(i?), 

k=i 

which occurs in our last definition, is valid if and only if ai, ct2, ■ ■ ■ , a m satisfy the following system 
of linear equations: 

m 

Y J ( X k) i ■ a k = i\ Xi (0<i<s-l). 

k=i 

From this we see that, for any given x G C s and m > s, it is possible to find X k ,ctk such that the 

above relation is valid. Consequently, \x\ is a well defined real number. 

In part 2 of the proof we show that \x\ is a norm for the vectors x £ C s . 

2. a) Let x,y G C s , and J> fc e AfcB = EAe wB = We have 

k z 

^ a fc e A * B + £ A e wB = (x + 

and therefore 

\x + y\ < J^lafclel^l +^|A|el w L 
k i 

This implies that 

\x + y\ < \x\ + |y|. 

b) For x G C s , A G C and any a fc , A fc G C with J2 a k e XkB = x(B), 

k 

we have ]T(Aa fe ) e AfcB = (Xx)(B), so that 
k 

\Xx\ < ^|Aa fc |e |Afc| = |A| -^|a fe |e |Afel . 

k k 

This implies that \\x\ < |A| • \x\. For A / we thus have \x\ = |A _1 (Ax)| < |A _1 | • |Ax|, and therefore 
|Ax| > |A| • \x\, which implies 

|Ax| = |A| • \x\. 

The last equality is also valid for A = 0. 

c) Suppose \x\ =0. For each e > there are a k , X k with 



^2a k e XkB = x(B) and ^|a fc |e |A,t| < e. 



k k 
For all such a k , X k and < i < s — 1 there follows 



J> fc (A fc )7*! < EKI^ ^ E|a fc |el A ^l<e. 

k k k 



42 



Consequently, \x{\ = 0, and therefore x = 0. 

d) We have proved that \x\ is a norm for x G C s , and we denote the corresponding 
matrix norm on C s,s by || • ||. 

3. In order to prove (i) we consider arbitrary x,y G C s and (eC with 

y = e^ B i. 

We have 

i n 

y i = S ^ — x i - j (i = 0, 1, . . . ,s - 1), 
j=o J ' 

and therefore also 

y(B) = e< B x(B). 

For any ctk, A& G C with ^ e XkB = x(B) there follows y(B) = ^ a*, e^ B e Afc ' B . In view of 

k k 

Theorem 3.3.1(c) we obtain 

y(B) = ^kS + ^ B . 

k 

Consequently 

\y\ < ^|« fc |e K+Afc| < e^^|a fc |el Afc l. 
k k 

It follows that \y\ < e'^'|x|. Hence 

||e CB || < e lcl for all ( 6 C. 

The matrix B thus satisfies condition III of Theorem 3.5.3 (with A = B, and with M = 1, 
W = {( : |C| < 1}). By virtue of that theorem, condition II is satisfied as well. This proves (i). 

4. Let v £ C s , with components vo = 1, Vi = (1 < i < s — 1). Let 1 < n < s — 1 and 
define the vector w G C s by 

u> = B n v. 

We shall show that 

M < 1, 

and 

M > n! (e/n) n . 

In view of the Theorems 5.3.1 and 2.1.2, the last two inequalities prove (ii). 

m 

The first of these inequalities follows from the fact that Yl a k e XkB = v(B) holds with 

k=i 

m = 1, oil = 1, Ai = 0. 

In order to prove the second inequality we consider arbitrary A& such that ak e XhB = 

k 

w(B). Since the components of w are given by 

Wi = (0 < i < s — 1, i / n), w n = 1, 
we have w(B) = B n , and there follows 



fc 



n! ^-^ 
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|A fc | n e-l Afe l n n e- n 

= i — i ' 

n! n! 

n!(e/n) n < ^|a fc |e |Afc| , 
k 

□ 

5.3.2 Improvements of (5.3.1) for the case of special norms 

According to the following theorem the stability estimate (5.3.1) can be improved substantially for 
the case of some important matrix norms. 

Theorem 5.3.3 Let s > 1, Q G C s ' s invertible, and p = 1,2 or oo. Let the norm \\ ■ \\ on C s,s he 
defined by \\A\\ = WQAQ' 1 ^ (for all A G C s ' s ). Then (5.1.2) with M 1 = 1 implies that 

\\B n \\ < M ( for n = 0,1,2,...), 

with M = 1 (if p = 1 or oo) or M = 2 (ifp = 2 ). 

Proof. Since the result for general invertible Q easily follows from the result for Q = I, it is sufficient 
to consider the latter case only. 

Let p = oo. Suppose B = (Pjk) satisfies (5.1.2) with || • || = || • ||oo, M\ = 1. Clearly (II) 
(Section 3.5) holds with W = {( : \(\ < 1}, M = 1. By Theorem 3.5.3 we have t^B, 1] C W. In 
view of Theorem 3.5.5 the matrix B satisfies a circle condition with respect to the unit disk, which 
means that ||-B||oo < 1- Therefore (5.1.1) holds with M$ = 1. 

For p = 1 the proof follows from the result for p = oo and the fact that \\A\\i = \\A*\\ 00 for 
all A G C s ' s . 

For p = 2 the proof runs as follows. Similarly as above it can be seen that (5.1.2) implies 
t 2 [B, 1] C W = {( : \C\ < 1}. In view of Theorem 3.5.4 we thus have 

{x*Bx : x G C s with x*x = 1} C W. 

The proof continues by applying Berger's inequality. This inequality reads 

r(A n ) < [r(A)] n for n= 1,2,3,... 

where A is any s x s matrix, and r(A) denotes the so-called numerical radius of A defined by 

r(A) = max {|ir*Ac| : x G C s with x*x = 1}. 

Since r(B) < 1, there follows 

r{B n ) < 1. 

We split B n into a sum B n = A\ + i^2 with Hermitian matrices A\ = ^(A + A*) and 
A 2 = ^ {A ~ A*)- By noting that for any Hermitian A the relation ||A|| 2 = r(A) is valid, we finally 
obtain 

\\B n \\ < Pill + p 2 || =r(A 1 ) + r(A 2 ) < 2r(B n ) < 2. □ 



where 



Consequently, 



and therefore \w\ > n\ (e/n) n . 
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5.4 The best stability estimates for fixed Mi > 1 

Theorem 5.2.1 shows that if the resolvent condition (5.1.2) is satisfied with a fixed Mi, then ||i? n || 
can grow at most linearly with n or s. Corollary 5.2.3 reveals that the corresponding upper bound 
is sharp — if we allow M 1 to be variable. 

For the special case Mi = 1, however, this linear growth with n or s is too pessimistic, as 
can be seen from the Theorems 5.3.1 and 5.3.3 in the previous section. 

For other fixed values Mi > 1, also, the question arises of whether the upper bound (5.2.1) 
can be improved. According to Theorem 5.3.2, a growth of \\B n \\ at the rate y/n or y/s can occur 
if arbitrary norms are considered. But, that theorem is not relevant to the important norms || • || p , 
with p = 1, 2, oo. In the following we present two results for these norms. The first result shows, 
for p = 1, oo, that growth at the rate n and s can occur, whereas the second result establishes, for 
p = 2, a growth which is almost at the rate y/\ogn and \/log s. 

Theorem 5.4.1 Let Mi > tt + 1, and p = 1 or p = oo. Let s be a given integer, with s > 1. Then 
there is an s x s matrix B such that (5.1.2) is satisfied with \\ ■ \\ = \\ ■ \\ p , and 

\\B n \\ p >n for n = 1,2, ... ,s. 

Theorem 5.4.2 Let Mi > tt + 1 be given. Then there exist a constant C > and matrices 
B s G C s,s for s = 2, 4, 6, . . ., such that all B s satisfy (5.1.2) with \\ ■ \\ = \\ ■ \\ 2 , and 

log(logs) 

Proof. By McCarthy & Schwartz (1965) it was shown that a constant 7 > and s x s matrices 
E s j (for all even positive s and j = 1, 2, . . . , s) exist with the following properties: 

s 

(5.4.1) !■:,.,/()■ E sJ E s , k = O (j + k), Y / E aJ =I; 

3=1 



(5.4.2) 



j odd 



> 7(logs) 1/2 / log logs; 



(5.4.3) B s = ^2e 2nij/s E SyJ satisfies (5.1.2) . 

j=i 

For even s we have 

<">r 2 ^ / 2^ /•,,. 

j = l j odd 

In view of (5.4.2) this implies 

||(£ s r /2 || 2 > -l + 2 7 (logs) 1 / 2 / log logs for s = 2,4,6,... . 
Since all (B s ) s / 2 / there exists a constant C with the property stated in the theorem. 
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5.5 Notes and remarks 

The Kreiss matrix theorem was originally published in Kreiss (1962). Subsequent elaborations and 
applications to so-called amplification matrices, originating from Fourier transformation, can be 
found e.g. in Richtmyer & Morton (1967), Strikwerda (1989). 

As already mentioned in Section 5.1, the Kreiss matrix theorem asserts, for the spectral norm, 
that the resolvent condition (5.1.2) implies power boundedness (5.1.1) with a stability constant M 
depending only on Mi and the dimension s. According to Tadmor (1981), the original proof by 
Kreiss (1962) yields an upper bound ||-B n || < M with 

M ~(M 1 ) sS , 

which is far from sharp. After successive improvements by various authors (cf. Morton (1964) and 
Miller & Strang (1966)), it was Tadmor (1981) who succeeded in proving a bound that is linear in 

s, 

\\B n \\ < 32evr- 1 sM 1 . 

LeVeque & Trefethen (1984) lowered this upper bound to 2esMi, and conjectured that the latter 
bound can be improved further to (5.2.9) (with || • || still standing for the spectral norm). Moreover, 
these authors showed by means of a counterexample that the factor e in (5.2.9) cannot be replaced 
by any smaller constant if the upper bound should be valid for arbitrary factors Mi in (5.1.2) and 
arbitrarily large dimensions s. 

Smith (1985) proved a result which, combined with the arguments of LeVeque & Trefethen 
(1984), leads to the bound ||-B n || < 7r _1 (7r + 2)esMi, which is an improvement over the upper bound 
2esMi but still weaker than (5.2.9). Eventually (5.2.9) was proved to be true (in Spijker (1991)). 
For an interesting historical survey see Wegert & Trefethen (1994). 

The proof in Section 5.2 of (5.2.1) has been taken from Dorsselaer, Kraaijevanger & Spijker 
(1993), and is partly based on arguments used earlier by Lenferink & Spijker (1991a,b), Lubich & 
Nevanlinna (1991), Reddy & Trefethen (1990). The crucial idea to use in the proof a relation of 
the form (5.2.5) and to bound the integral f r |-R'(C)| \d(\ in terms of max \R(()\ was used earlier 

by LeVeque & Trefethen (1984). 

Corollary 5.2.3 was proved by LeVeque & Trefethen (1984) for the spectral norm. The proof 
of Theorem 5.2.2 has been taken from Dorsselaer, Kraaijevanger & Spijker (1993). 

The proof of (5.3.1) in Section 5.3.1 is essentially based on ideas taken from Bonsall & 
Duncan (1980)(see also Bonsall & Duncan (1971)). Another proof can be given along the lines of 
Lubich & Nevanlinna (1991) (Theorem 2.1). 

Our proof of Theorem 5.3.2 strongly relies on ideas taken from Crabb (1970). 

The value Mq = 2 (for p = 2) in Theorem 5.3.3 has been known already for some time 
and was stated e.g. in Reddy & Trefethen (1992). For the inequality of Berger, used in the proof 
of Theorem 5.3.3, we refer to Bonsall & Duncan (1980), Horn k Johnson (1990), Pearcy (1966), 
Richtmyer & Morton (1967, p. 89). 

Theorem 5.4.1 follows easily from a clever counterexample presented in Kraaijevanger (1994). 

Theorem 5.4.2 and its proof have been taken from Dorsselaer, Kraaijevanger & Spijker 

(1993). 
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(6.1.1) 



6 Stability estimates under resolvent conditions on hA 
6.1 Linear stability analysis and stability regions 

Consider an initial value problem for a system of s ordinary differential equations of the form 

(U'(t) = AU{t) + r{t) (t>0), 
\U(0) = n . 

Here A is a given constant s x s matrix, and no, r(t) are given vectors in C s . The vector U (t) G C s 
is unknown for t > 0. 

In this section we analyse the stability of numerical processes for approximating U (t) . This 
analysis will be relevant also to classes of numerical processes for solving partial differential equa- 
tions. 

To elucidate this relevance, we assume an initial-boundary value problem to be given for a 
linear partial differential equation with variable coefficients (which depend on the space variable x 
but not on the time variable t). Applying the method of semi-discretization, where discretization 
is applied to the space variable x only, one can arrive at an initial value problem for a large 
system of the form (6.1.1). In this case the matrix A, the inhomogeneous term r(t), and the 
vector no are determined by the original initial-boundary value problem and by the process of 
semi-discretization. The solution U(t) to (6.1.1) then provides an approximation to the solution 
of the original initial-boundary value problem. For examples we refer to the Sections 1.2, 3.6, 4.2, 
where semi-discretization by finite- difference methods is dealt with . 

Many step-by-step methods for the numerical solution of ordinary differential equations, like 
Runge-Kutta methods or Rosenbrock methods reduce — when applied to (6.1.1) — to processes 
of the form 

(6.1.2) u n = (p{hA)u n -\ + r n for n= 1,2,3,... . 

Here (p(z) = P(z)/Q{z) is a rational function, depending only on the underlying step-by-step 
method, and P (z), Q(z) are polynomials, without common zeros, such that <^(0) = <p'(0) = 1. 
Further, h = At > is the stepsize, and ip(hA) = P(/iA)[<5(/iA)]" 1 (when Q(hX) / for all 
A G <t[A], cf. Example 3.3.2). The vectors r n G C s are related to r(t), and u n ~ U(nh) are 
computed successively from (6.1.2). It is worth noting that many numerical processes in partial 
differential equations which are not constructed with the process of semi-discretization in mind are 
still of the form (6.1.2), and can a posteriori be conceived as relying on semi-discretization. 

An example of (6.1.2) is provided by the, fully discrete, numerical process constructed in 
Section 4.2 for the solution of (4.2.2). From (4.2.3) we see that this process is of the form (6.1.2) 
with 

<p(z) = (i + (i-0)z) (l-ez)- 1 

and the tridiagonal matrix A = (ccjfc) given by (4.2.4). 

Since (6.1.2) is a special case of (4.2.1) the stability analysis of (6.1.2) amounts, as explained 
at the end of Section 4.2, to investigating the growth of the matrices B n specified by 

(6.1.3) B = ip(hA). 

In this analysis it is useful to introduce the stability region S, defined by 

(6.1.4) S = {z: z G C with Q(z) / and \(p(z)\ < 1}. 

The following theorem is a variant to Theorem 4.3.1. It is easier to apply than the latter, 
since it gives a stability criterion in terms of <p(z) and hA rather than in terms of the (more 
complicated) matrix B itself. 
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Theorem 6.1.1 (The eigenvalue criterion) Let \ \ ■ \ \ be an arbitrary induced norm on C s ' s , and 
B = ip(hA). Then there is an M such that (5.1.1) holds, if and only if 



(6.1.5) 



all eigenvalues of hA belong to S; and whenever Jk is a Jordan block 
(with order Sk > 1) of hA corresponding to an eigenvalue h\k £ dS, 
, then the derivatives tp^ (h\k) vanish for j = 1, 2, . . . , Sk — 1. 



Proof. Applying Theorem 3.3.6 (with /, A replaced by tp, hA) we see that 

r 

B n = <p(hA) n = £(T fc ) n , 



fc=i 



where 



T fc = <p(h\ k )P k + Qk, Qk=Yl * , k) (RkY, A fc G Sfe > 1. 

J=l J ' 

There is an Mq with (5.1.1) if and only if, for each fc, the powers (Tk) n remain bounded for n — > oo. 
Let be given. From the above expressions for T^, we see that 

(T fc )" = [<^(/iA fc )P fe + Q k ] n = ( n Mh\ k ) n - p Ql P k , 

where nik = min(n, — 1). Moreover, 

Qk = O if and only if <p^\h\k) = for 1 < j < Sk — 1. 
The proof of the theorem is completed by noting that the following four implications are 

valid. 

(i) \ip(h\k)\ < 1 =>■ (Tk) n remains bounded for n — > oo; 

(ii) |(^(/iAfc)| = 1 and Qk = O (Tfc) n remains bounded for n — > oo; 

(iii) \ip(h\k)\ = 1 and Qk ^ O (Tfc) n does not stay bounded for n — > oo; 

(iv) |i^(/iAfc)| > 1 (Tfc) n does not stay bounded for n — > oo. □ 

Many functions (^(z) of practical interest have non- vanishing derivatives <p'(z) on the whole 
of dS. In this case (6.1.5) simply reduces to a[hA] C S and the condition that the Jordan blocks 
Jk of hA have order 1 whenever they correspond to an eigenvalue h\k € 9»S. 

Clearly, condition (6.1.5) always implies that 



(6.1.6) <r[/iA] C 5. 

Moreover, from Theorem 3.1.4 we see that (6.1.6) is equivalent to (6.1.5) if the matrix A is normal. 
The analogue of Theorem 4.3.2, in the situation where B = ip(hA), is as follows. 

Theorem 6.1.2 ( The eigenvalue criterion for normal matrices) Let \ \ • \ \ denote the spectral norm, 
and assume B = ip(hA) where A is normal. Then the following three statements are equivalent to 
each other. 

(i) There is an M such that (5.1.1) holds. 

(ii) \\B n \\ < 1 for n = 1,2,3,.... 
(Hi) a[hA] c S. 
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Proof. Suppose (5.1.1) holds for some constant M . Then, in view of Theorem 6.1.1, condition 

(6.1.6) is fulfilled. 

Suppose (6.1.6) holds. By Theorem 3.1.4, all Jordan blocks of hA have an order 1, and 
Theorem 3.3.6 thus yields the representation 

r 

B = p(hA) = Y,v{h\ k )Pk =Tdiag [^(/iAi), <p(h\ 2 ), ■ ■ ■ , <p(h\ s )]T-\ 

k=l 

Consequently ||-B n || < ||T|| • = 1 as in the proof of Theorem 4.3.2. □ 

In general (6.1.5) has similar advantages and disadvantages as the eigenvalue condition 
(4.3.1). It is relatively simple to verify, and reliable in the situation of normal matrices, for which 
it reduces to (6.1.6). But, it is unreliable for (families of) matrices that are not normal. 

In the rest of this chapter we adapt (6.1.5) to conditions on hA that reliably predict 
stability — also for nonnormal matrices and norms || • || on C s ' s different from the spectral norm. 
An advantage of these conditions on hA over a resolvent condition on B = <p(hA) (as dealt with in 
Chapter 5) lies in the circumstance that, in general, hA has simpler a structure than B, and that 
knowledge available about S can be exploited. 

The conditions on hA we shall deal with, are all of the following general form, 

(6.1.7) (zl -hA) is regular and || (zl - hA)' 1 1| < K ■ d(z, V)~ l 
for all complex numbers z G C \ V. 

Here K is a constant, V a closed subset of C, || • || denotes an arbitrary induced norm on C s ' s , and 
d(z, V) is the distance from z to V. 

In case V is bounded, we easily see from (6.1.7), by letting z — > oo, that 

(6.1.8a) K > 1. 

In the following we shall always assume that this inequality is fulfilled. 
We further note that if 

(6.1.86) V is a closed subset of S, 

condition (6.1.7) implies (6.1.5). We shall focus on the situation where (6.1.8b) is fulfilled. 

In the following sections, condition (6.1.7) will be seen to imply stability estimates of the 

form 

(6.1.9) \\cp(hA) n \\ < K -<S>(n,s) for n = 1,2,3, ... , 

where the function 3> only depends on ip and V (and not on h,A,K or || • ||). 

6.2 Stability estimates which grow linearly with n and s 

6.2.1. Arbitrary subsets V of the stability region 

The following general theorem is valid. 

Theorem 6.2.1 There is a constant 7 such that the stability estimate (6.1.9) holds, with 

<&(n, s) = 7 • min{n, s}, 

whenever K, V satisfy (6.1.8) and hA G C s ' s satisfies condition (6.1.7). Here 7 only depends on 
the rational function tp(z) (and not on V, n > 1, s > 1, hA G C s ' s , || • || or K). 
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This theorem will be proved below for the special case where 

(6.2.1a) ip'(z) / on dS, 
and 

(6.2.1b) S is bounded. 

The proof will consist of two steps. First, it will be shown that B, given by (6.1.3), satisfies a 
resolvent condition (5.1.2) (with the same norm || • || as appearing in (6.1.7), (6.1.9)). A subsequent 
application of Theorem 5.2.1 (part (b)) will lead to (6.1.9) with <3?(n,s) as specified in Theorem 
6.2.1. 

In proving that B satisfies (5.1.2) we shall make use of two lemmas, the first of which reads 
as follows. 

Lemma 6.2.2 Let p > 1, M > 0, B € C s,s and || • || an induced norm on C s ' s . Assume 

(6.2.2a) ((I - B) is regular for all ( with \(\ > 1, and 

(6.2.2b) - -B) -1 || < M ■ (|C| - l)" 1 for all ( with K |C| < P- 

Then (5.1.2) is satisfied with 

Mi = ^- M. 

Proof. For any C with \(\ > p we can write ((I — B)~ l as a Dunford- Taylor integral, 

(CI - B)- 1 = ^-J r f(z)(zl - B)- l dz, 

where f(z) = (( — z) -1 and T is any positively oriented circle r[0, a] with 1 < a < p. Consequently, 
for such C, we obtain 

m-B^w < (27T)- 1 / \c - z \-i • m • (\ z \ - i)-'\dz\ < M ; CT 

Choosing a = ^fp we obtain, for \C\ > p, 

(KI _ 1)ll(a _ Brl|l ^ (l+ ^l),^. M . 

Since Mi > M the proof is complete. □ 

For a > 0, (3 > we shall use, in the following, the notations 

D a ={(: CGCwithK|C|<l + a}, 
Sp = {z : zeC with < d(z, S) < (3}. 

For given (eCwe shall deal with the rational function 

(6.2.3a) ^(z) = [(-<p(z)}- 1 . 
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It is easily verified that the order of tp^(z) does not exceed the order, say r, of <p(z). Therefore the 
number of different poles of V"c( z )> denoted by q((), satisfies 

(6.2.36) q(Q < r. 

We denote the poles of ip((z), an d their corresponding orders, by 

(6.2.3c) Zj (0 and kj(Q ( for j = 1, 2, . . . , q(Q), 

respectively. Finally, we use the notation 

(6.2.3d) ip(oo) = lim ip(z). 

z— >oo 

Our second lemma only concerns (f(z) (and not the matrix B given by (6.1.3)): 

Lemma 6.2.3 Let tp(z) satisfy (6.2.1a,b). Then there are positive a, (3 with the following proper- 
ties. 

(i) \ip(oo)\ > I + a; 

(ii) ip(z) is regular, and if'(z) / at all z G Sp; 

(iii) if C G D a , then Zj (() G Sp (for j = l,2,..., g(Q); 

-1 -1 

(iv) if C G £>„, then V<(*) = [C - ^M]" 1 - £ [^(C))] (* - *;(0) ; 

.7 = 1 

(v) if C G £> a , then |C| - 1 < Mi ■ d(z 3 -(C), 5) (for j = 1, 2, • • • , g(C)j, 
with n.\ = max{|(/2'(z)| : z G cl (S^)} < oo. 

Proof. 1. Since S 1 is bounded, <p(oo) satisfies 1 < |<^(oo)| < oo. In case |y(oo)| > 1, it is clear that 
(i) holds for all a > sufficiently small. 

Suppose y?(oo) = rj with \rj\ = 1. Then we can write 

P(-z) = a r z r + • • • + aiz + ao, Q{z) = f3 r z r + • • • + faz + (3 , 

where a r = n- (3 r ^ 0. Therefore, for w G C with |«;| > sufficiently small, we have a representation 
for ip(l/w) of the form 

1 ~t~ ^j\w ~\~ ' ' ' ~\~ ^y r iv r 

ip(l/w) = rj ■ f(w) with f(w) 



1 + 8\w + • • • + 5 r w r 
Since 

f{w) = 1 + ei ■w k + 0{w 2 ) ( for w -►()), 

with ei / 0, > 1, it follows that there exist w, arbitrarily close to 0, with \f(w)\ < 1. This means 
that ip(l/w) assumes values with modulus less than 1 for w arbitrarily close to 0. This contradicts 
(6.2.1b). 

We conclude that (i) holds for all a > that are sufficiently small. 

2. Denote the set of all poles of (p(z) and zeros of <p'(z) by T. By (6.1.4), (6.2.1a) the sets 
dS and T have no points in common. Since dS is closed and T finite, there exists a (3 > such 
that all points in T have a distance to dS greater than (3. 

Let z G Sp. Since d(z,dS) = d(z,S) it follows that < d(z, dS) < (3. Consequently z does 
not belong to T. Part (ii) has been proved. 

3. We define 

a = inf {\<p(z)\ : z G C with d(z, S) > /?}. 
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First we assume a > 1. We choose a satisfying (i) with < a < a — 1. Let £ be given, with 
C €. D a . We then have 1 < \C\ < o. The definition of a and (6.1.4) imply that all z with (p(z) = £ 
must belong to Sp. Hence Zj(Q 6 for j = 1, 2, . . . , q. 

Next we assume a < 1, which will be shown to lead to a contradiction. The definition of a 
implies that there is a sequence of complex numbers j/i, J/2, J/3) ■ ■ ■ with 

d(y k ,S)>P and |^(y fe )|<l + l/fc ( for fc = 1,2,3, .. .)■ 

From (i) we see that positive numbers p and e exist such that |¥>(z)| > 1 + e for all complex 
z not belonging to the disk D[0,p]. Hence 

yfe G D[0, p] for all k that are sufficiently large . 

Since D[0,p] is compact, there is a convergent subsequence of {?/&}, say {yk(j)} with limit yoo. 

Since d(y fc (j), S 1 ) > /3, we have d(y 00 ,S') > /3. But, since <p(yk(j)) tends to <f(yoo) (for 
j — ► oo), we also have |<p(yoo)| < 1- The last inequality means that d(y 00 ,S') = 0, and we have a 
contradiction. 

4. Let £ e D Q . In view of (ii), (iii) we see that at all poles Zj(C) we have (p(zj(Q) / oo and 
<p'{zj(Q) / 0- From (6.2.3a) it follows that the principal part of 4>((z) at Zj(() equals 

Subtracting from tpci z ) an its principal parts we obtain a function that we denote by w$(z). 
In view of Theorem 2.1.1 the function oj^(z) is a polynomial in the variable z. Since ip^(z) tends 
to [£ — 99(00)] _1 when z — > 00, we can conclude that 

lim loAz) = \C — <p(oo)] 1 . 

Clearly, uj^(z) is a polynomial of degree zero, and (iv) thus holds. 

5. Let £ G D a and y G 96*. We have 

ICI-i < 1^(0)1-1^)1 < Mzj(0) -<p(y)\- 

Denote the straight line segment connecting Zj(() to y 6 95 with |zj(C) — y| = d(£j (C), S 1 ) by L. 
Since L is contained in cl (Sp) we have 

b(^-(0) - <p(v)\ < m-\zj(Q-y\ = m-d(zj(C),s), 

with \x\ as in statement (v). This completes the proof of the lemma. □ 

We now complete the proof of Theorem 6.2.1, making use of the Lemmas 6.2.2 and 6.2.3. 

Proof of Theorem 6.2.1 for the case (6.2.1a,b). Assume (6.2.1), (6.1.7), (6.1.8). The matrix 
B = ip(hA) exists (in the sense of Example 3.3.2) since o^/lA] C V C S and the denominator Q(z) 
of ip(z) does not vanish on S. Moreover ((I — B) is regular for all ( with \(\ > 1. This follows from 
the fact that a[B] = <p(a[hA]) C {C : |CI < 1} (by Theorem 3.3.1 (d) and (6.1.4)). 

Choose a, (3 according to Lemma 6.2.3. Let Q £ D a . In view of Lemma 6.2.3 (iv) and 
Example 3.3.2 we have 

9(0 

(CI-B)- 1 = [C-^oo)]- 1 /-^^^^))]- 1 ^-^^)/)- 1 . 
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By taking norms we obtain 

1 i(0 

llk'-^'ll * l^l-d + aJ -W-T'^gllfeK' 7 -^"'" 

where 

lUo = min{|^'(2;)| : z£ cl (5/?)} 

is positive by Lemma 6.2.3 (ii) and (6.2.1a). 

The assumptions (6.1.7), (6.1.8b), in combination with Lemma 6.2.3 (v), yield 

11(^(0/ - hAy'W < K ■ d( Zj (Q, V)- 1 < K ■ diz.iC)^)- 1 < Kp 1 • (|C| - l)" 1 . 
In view of (6.2.3b) we thus obtain 

Ua-By'w < m.(ici-i)- 1 , 



with 

M = 

\(p(oo)\ - (1 + a) fi 



M = , . r + rK- IU 



Clearly, the assumptions of Lemma 6.2.2 are fulfilled with p = 1 + a. Therefore, our matrix 
B satisfies (5.1.2) with 

M 1 = V l±^ + ) [( |y(oo)| -I- a )~ l a + w" 1 -^. 
Vl + a — 1 

In view of (6.1.8a), the matrix B satisfies (5.1.2) with 

Mi = 71 • K, 

where 71 only depends on the rational function (p(z). 
An application of Theorem 5.2.1 (b) yields 

\\B n \\ < (l + l/n) n min{n+l,s} 7l ^ < 2e 7l • K ■ min{n, a}. 

Theorem 6.2.1 has thus been proved, under the assumption (6.2.1a,b), with 7 = 2e7i. □ 

6.2.2 An example in which V is a disk and S is bounded 

Let || • || denote an arbitrary induced norm on C s,s , and suppose A G C s ' s satisfies the circle 
condition 

(6.2.4) \\A-jI\\<p. 

In applications of Theorem 6.2.1 the following lemma will be helpful. 

Lemma 6.2.4 Assume (6.2.4). Then, for any h > 0, the matrix hA satisfies the resolvent condition 
(6.1.7) with 

V = D[h-/,hp] and K=l. 
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Proof. From (6.2.4) it follows that 

\\hA-hjI\\ < hp. 

Denoting the 1-numerical range of hA by r[hA] (see Definition 3.5.1) we thus have 

r[hA] c D[h-y, hp]. 

An application of Theorem 3.5.3 (with hA replacing A, and with M = 1, W = D[hj, hp]) completes 
the proof. □ 

We shall illustrate Theorem 6.2.1 in the numerical solution of the diffusion-convection- 
reaction problem (4.2.2). We define the tridiagonal s x s matrix A by (4.2.4). With regard to 
the coefficients a(x), b(x), c(x) we make the assumption (3.6.1). 

In Section 3.6 we showed that the matrix A satisfies a circle condition (6.2.4) with || • || = || • || oo 

and 

7 = ~2' p= 2 : a = 4 ° Hoc + 20 |b|oo + |c|oo- 
By virtue of Lemma 6.2.4 the matrix hA satisfies the condition (6.1.7) with 

(6.2.5) II • II = II • Hoc, K=l and 

2h h h 

V = D[-r ,r ], where r = ^-|a|oo + ^|o|oo + ^\c\oo- 

We consider the, fully discrete, process of Section 4.2 with a matrix B given by (4.2.3), 
(4.2.4). This amounts to (6.1.2) with ip(z) = (1 + (1 - 9)z) (1 - 9z)~ x . We consider a fixed 9 with 
< 9 < \. The corresponding stability region S = S$ is given by the formula 

S = D[-r,r] with r = — ^— for < 9 < 1/2. 

1 — 29 

Clearly, all of the assumptions of Theorem 6.2.1 will be fulfilled as soon as r < r, i.e. 

2h h h 1 

(6-2.6) ^"l°loo + i:\b\oo + ^Icloo < Y^29' 

Under this condition we thus have stability, (6.1.8) with $ as in Theorem 6.2.1. 
As a numerical illustration we consider problem (4.2.2) with 

a(x) = 1, b(x) = -10 3 , c(z) = 0, 

and the corresponding numerical process with 

« = * = ! = !. 

4' s 20 

In this situation the sufficient condition for stability (6.2.6) amounts to the stepsize restriction 

h<h = ~ 9.62 x 10" 5 . 

~ u 10 400 

Computer experiments show that with a stepsize h = 9 x 10~ 5 one has 

\\<p{hA) n \\ < 2.4 for n = 0,1,2,... . 
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The last inequality expresses a stable behaviour, which is amply in agreement with Theorem 6.2.1. 

It is instructive to compare the above stepsize restriction based on Theorem 6.2.1 with a 
naive use of the eigenvalue condition as formulated in Theorem 6.1.1. In fact, for a(x), b(x), c(x) 
and 5 as specified above all eigenvalues A of A can be shown to be different from each other and to 
satisfy 

-26 514 < A < 0. 

For 9 = 1/4 we have S = D[—2, 2], so that the eigenvalue condition (6.1.5) is satisfied for all h > 
with 

" £ *■ = 2<rb = 15 1 x 10 " 

Computer experiments show that with a stepsize h = 15 x 10~ 5 one has 

max ||(^(/iA) n || ~ 2.7 x 10 12 . 

n>0 

Hence from a practical point of view there is mere instability with this h, although h < hi and 
(6.1.5) is fulfilled. 

6.2.3 An example in which V = C_ C S 

In the example of Section 6.2.2 we have seen that, in the situation (6.2.5), condition (6.1.8b) boils 
down to a stepsize restriction of the form h < ho-, with finite ho- This is related to the fact that 
the stability region S, dealt with in the example of Section 6.2.2, is bounded. In cases where S is 
unbounded it may be possible to establish stability estimates which are valid for all h > 0. 

In order to illustrate this point, we consider the tridiagonal matrix A given by (4.2.4) and 
assume (3.6.1) to be fulfilled. In Section 6.2.2 we have seen that hA satisfies the resolvent condition 
(6.1.7) with the specifications (6.2.5). Since D[—r ,r ] C C_ it follows that hA also satisfies 
condition (6.1.7) with 

|| • || = || • Hoo, K = 1 and V = C_ . 

We again consider the fully discrete process of Section 4.2 with B specified by (4.2.3), (4.2.4). 
But, now we consider a fixed 9 with \ < < 1. This amounts to (6.1.2) with a rational function 
ip(z) the stability region S of which satisfies 

C_ C S. 

An application of Theorem 6.2.1 shows that the process of Section 4.2, with \ < 9 < 1, is 
stable in the sense that 

||-B"1|oo < 7-min{n, s} for all n > 1, s > 1. 
Here 7 only depends on 9, and the result is valid for any h > 0. 

6.3 Stability estimates which grow slower than linearly with n 

Under conditions on V which are (slightly) stronger than (6.1.8b) variants to Theorem 6.2.1 exist 
in which 

3>(n, s) = 7-min{n Q ,s} 

with a < 1. Below we formulate such a variant with a = 0. 
We consider the situation where 

(6.3.1a) < 9 < 6' < vr/2, 

(6.3.16) W{9') C S, 

(6.3.1c) V = W{9). 

The following interesting theorem is valid. 
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Theorem 6.3.1 Assume (6.3.1). Then there is a constant 7 such that (6.1.7) implies (6.1.9) with 
s) = 7. Here 7 only depends on ip(z) and V (i.e. on 9). 



6.4 Resolvent conditions and the M-numerical range of hA 

6.4.1 The construction of a set V as in the resolvent condition (6.1.7) 
by using M-numerical ranges 

Let || • || be an arbitrary induced norm on C s ' s , and K > 1. Suppose A G C s,s , h > and 

W = r[hA,K] 
(cf. Section 3.5). From Theorem 3.5.3 we see that 

j zI-hA is regular and \\(zl - hA)~ l \\ < K • [d(z, VF)]" 1 
\ for all z W. 

Therefore, we can make the following two observations. 
(I) If V is any set with 

r[hA,K\ C V C C, 

then hA satisfies the resolvent condition (6.1.7). 
(II) In order to construct a set V as in observation (I) we only have to determine a finite number 
of pairs 7^, pj such that 

\\(hA-~f j I) k \\ < K( Pj ) k for k = 1,2,3,... . 

In view of Definition 3.5.1 the set 

V = f\D[y j ,p j ] 

j 

is as required. 

6.4.2 An illustration in the numerical solution of the pure diffusion equation 

We consider an initial value problem of the form (6.1.1) originating from a semi-discretization of 
the pure diffusion problem(1.1.4). We assume the semi-discretization to be based on the finite 
difference formula (1.2.2a). The s x s matrix A = {ctjk) in (6.1.1) thus equals the tridiagonal 
matrix A displayed in Subsection 1.2.3, i.e. 



(6.4.1) 



/ a jj = - 2 /$ 2 i a jk = V<5 2 for \j - k\ = 1, a jk = for \j - k\ > 1, 
\ and 5 = (s + l)" 1 . 

Without proof we state the following result on the matrix A. 



Lemma 6.4.1 Let a be given with < a < tt/2. Then there exist constants K = K a > 1 and 
R = R a > such that, for all s > 1, the s x s matrix A = (ajk) given by (6.4.1) satisfies 

IKl + e^-^R-WA^W < K for fc = 1,2,3,... . 
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From this lemma we easily see that, for k = 1, 2, 3, . . . , 



||(^-7i/) fc |L < K a (pi) k with 7l = -Ai ?Q e i ^-f), Pl = | 7 i|, 

and 



5 2 

IKhA-^I^W^ < K a (p 2 ) k with 72 = ( 7l )*, p 2 = |72| 



In view of the above observations (I), (II) we deduce that, for < a < 7r/2, the matrix hA 
satisfies (6.1.7) with 

K = K a and V = V a = L>[ 7l , Pl ] f| D[ 72 ,p 2 ]. 
We note that the set V a satisfies 

V a C W(a). 

This inclusion implies the important fact that for any a £ (0, vr/2) there is a constant .ftT = X Q 
such that the matrix given by (6.4.1) satisfies (6.1.7) with V = W(a). Here K a is independent of 
s > 1 and h > 0. Combining this fact with Theorem 6.3.1 we arrive at the following theorem. 

Theorem 6.4.2 Let (p(z) be such that, for some (5 G (0,7r/2), the wedge W{(3) is contained in the 
stability region S. Let the s x s matrix A = (etjk) be defined by (6.4.1). Then there is a constant 
7 such that \\(p(hA) n \\ 00 < 7. Here 7 only depends on ip(z) (and not on n > 1, s > 1 or h > 0). 

6.4.3 An illustration in the numerical solution of a diffusion-convection- 
reaction problem 

In the following we illustrate the relevance of observation (I) of Subsection 6.4.1 with K = 1 and 
We consider the initial-boundary value problem 



u(0,t) = go(t), u(l,t) = gi(t), u(x, 0) = f(x), where 0<z<l and t > 0. 



(6.4.2) —^u(x, f) = — u(x, t) - 200— u(x, t) - 137 000 • x ■ u(x, t), 



Here go(t), 9i(t), f(x) are given functions, and u(x,t) is unknown. We apply the method of 
semi-discretization using the finite difference formulas (1.2.1c), (1.2.2a). In this way (6.4.2) is 
transformed into a system of ordinary differential equations of the form (6.1.1) with a tridiagonal 
s x s matrix A = (ajk) for which 



(6.4.3) 



ctjj-! = <T 2 + 100 5' 1 for 2 < j < s, 
a jtj = -2<T 2 - 137 000 ■ j ■ 5 for 1 < j < s, 
, atj tj+1 = <T 2 - 100 5' 1 for 1 < j < s - 1, 



with 5= (a + 1) -1 . 

With this matrix A we first consider the numerical process (6.1.2) where 

<p(z) = 1 + z + 0.5 z 2 + 0.0625 z 3 . 



The stability region S corresponding to this function ip(z) contains the real interval [—6,0]. 
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Let s = 99, so that 8 = 10 -2 . Now all eigenvalues A of A are different from each other, and 
real with —157 000 < A < —20 000. Therefore the eigenvalue condition (6.1.5) is fulfilled whenever 
the stepsize ft > satisfies 

ft < ~ 3.82 x 10 -5 . 

~ 157 000 

We choose ft = 3.4 x 10~ 5 so that (6.1.5) is amply fulfilled. But, straightforward numerical experi- 
ments show that 

(6.4.4) max \\ip(hA) n ^ > 3 x 10 11 (for s = 99, ft = 3.4 x 10" 5 ). 

n>l 

This inequality once more illustrates the fact that Theorem 6.1.1 can be an unreliable guide to 
stability. 

Next we consider, with the same matrix A, the numerical process (6.1.2) where 

ip(z) = 1 + z + 0.5 z 2 + 0.0645 z 3 . 
For a > 0, p > we introduce the subset V(a, p) of the complex plane defined by 

V(a, p) = {z : z = x + y with x£K, y G C, —a — p < x < — p, \y\ < p}. 
With (To = 4.67 and po = 0.68 the set V(a ,p ) is contained in the stability region of <p(z), i.e. 

V(a ,p ) C S. 

Using Theorem 3.5.4 we easily see from the expressions (6.4.3) that, for any 5 > with 
5 < 10~ 2 , the 1-numerical range of A satisfies 

[A, 1] c V (137 000 , 2<T 2 ). 

In view of Theorem 3.5.2 (part (iii)) we thus arrive at 

Too[hA,l] C V (137000/i, 2h5~ 2 ). 

so that observation (I) of Subsection 6.4.1 applies to the situation at hand. 
By virtue of Theorem 6.2.1 we can conclude that 

(6.4.5) ||^(/iA) n || < 7-min{n,s} for all n > 1, s > 99 and h with 

0< ^- mi H 137000' Y' P0 \- 
Here 7 is a constant independent of n, s, h. 

Let s = 99. Since now 5 = 10~ 2 and therefore 

minf— ^— ,-.p ) = 3.4x 10" 5 , 
1137 000 2 J 



the stepsize restriction in (6.4.5) is fulfilled when 

ft = 3.4 x 10" 5 . 

With this stepsize we thus expect a mild error propagation in accordance with the stability estimate 
in (6.4.5). Straightforward numerical experiments show that actually 

(6.4.6) H^M^Iloo < 3/2 ( for all n > 1, s = 99, ft = 3.4 x 10 -5 ) 

which is in agreement with (6.4.5). 

The striking difference between (6.4.4) and (6.4.6) is closely related to the fact that the 
stability region 5 corresponding to (6.4.4) does not enclose the set T^hA, 1]. In fact, this S is so 
'small' that it does not contain any complex number z with Re(z) = —4 and lm(z) / 0. 
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6.5 Resolvent conditions and the e-pseudospectra of hA 

In the above we have seen that the resolvent condition (6.1.7) is a handy tool for deriving stability 
estimates of the general form (6.1.9). Further, we have seen that circle conditions and M-numerical 
ranges can provide sufficient conditions in order that (6.1.7) is fulfilled. But up to now an easy 
interpretation of (6.1.7) has been missing. 

We formulate a theorem which shows that (6.1.7) can nicely be interpreted in terms of the 
e-pseudospectra of the matrix hA. 

Theorem 6.5.1 Let V be a closed subset of C. Then the resolvent condition (6.1.7) is equivalent 
to the requirement that, for each e > 0, the set u £ [hA] is contained in 

{z : z eC with d(z, V) < K ■ e} = {z : z = x + y with x € V, \y\ < Ke}. 

This theorem can be proved in a straightforward way by using the fact that, according to 
Theorem 3.2.3, the statements (i) and (v) of Section 3.2 (with B = hA) are equivalent. 

We note that the concept of e-pseudospectra can also be used, in principle, to determine 
numerically regions V and constants K such that (6.1.7) holds. In order to explain how this may 
be done we write B = hA, choose a fixed e > and denote the boundary of a £ [B] by T £ . The set 

(6.5.1) V = o £ [B] 

can be determined numerically, e.g. by checking for a large set of complex numbers A whether (v) 
(Section 3.2) is satisfied. A corresponding constant K can be computed from the formula 

(6.5.2) K = \T £ \ (27re) _1 , 

where |r e | denotes the length of r e . 

In order to establish (6.5.2) we note that for z V we have 



{zI ~ B) ~ 1 = hJ T (^- A )" 1 ( A/ -^)" 1 



and therefore 

IKzI-B)" 1 !! < ^ max ((z-A)" 1 ^" 1 = K d (z, V) -1 . 

It is clear that both V and K depend on e. Therefore, it may pay to evaluate (6.5.1) and 
(6.5.2) for various values of e > 0. 



6.6 Notes and remarks 

We note that problems of the form (6.1.1) can originate form partial differential equations in cases 
where the process of semi-discretization does not rely on the introduction of finite differences, but 
instead on e.g. finite volumes, finite elements or (pseudo) spectral approximations (cf. Section 1.3). 

Step-by-step methods for the numerical solution of initial value problems in ordinary dif- 
ferential equations are discussed e.g. in Butcher (1987), Hairer & Wanner (1996). Runge-Kutta 
methods as well as Rosenbrock methods are discussed extensively in these works. 

Condition (6.1.5) was stated earlier in Dorsselaer, Kraaijevanger & Spijker (1993). 
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The proof of Theorem 6.2.1 for the special case (6.2.1), as presented in Section 6.2.1, is 
essentially due to Reddy & Trefethen (1992). The value for Mi, given in Lemma 6.2.2, has been 
taken from Spijker (1996). 

Lubich Sz Nevanlinna (1991) were the first to prove a version of Theorem 6.2.1 in which 
condition (6.2.1b) is not required; they proved (6.1.9) with 3>(n, s) = 7 • min{n, s} for the case 
where (6.1.7) holds with V = C_ C S. A proof of Theorem 6.2.1, for the general case, can be 
found in Spijker & Straetemans (1996b). 

The (numerical) example in Section 6.2.2, with tp(z) = (1 + (1 — 6)z) (1 — 9z)~ l , 9 = 1/4, 
has been taken from Kraaijevanger, Lenferink & Spijker (1987). 

Interesting stability estimates were derived, independently of each other, by Palencia (1993, 
1995) and Crouzeix, Larsson, Piskarev k, Thomee (1993). These estimates, adapted to fit in 
our terminology, yield Theorem 6.3.1. For further variants to Theorem 6.2.1 with 3>(n, s) = 
7 • min{n Q , s}, < a < 1, we refer to Spijker & Straetemans (1996b). 

Part of the material in the Sections 6.4.1, 6.4.2 has been taken from Lenferink & Spijker 
(1990), Spijker (1993). The proof of Lemma 6.4.1, as given in the first of these two references, 
relies on stability estimates of Thomee (1965). The example in Section 6.4.3 is due to Lenferink & 
Spijker (1991b). 

Material, closely related to Section 6.5, can be found in Dorsselaer, Kraaijevanger & Spijker 
(1993). Theorem 6.5.1 was formulated earlier in Reddy & Trefethen (1990, 1992). 
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